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ABSTRACT 


The  initiation  of  graupel  growth,  was  simulated 
by  approximating  a  hexagonal  ice  crystal  to  be  a  very 
thin  oblate  spheroid.  The  flow  of  air  passing  oblate 
spheroids  of  ice  with  Reynolds  numbers  of  10  and  20  was 
calculates  by  solving  the  steady  state  Navier-Stokes 
equation  in  oblate  spheroidal  co-ordinates  for  spheroids 
with  aspect  ratio  of  0.05.  The  flow  field  of  air  around 
the  water  drop  was  evaluated  by  two  methods:  solving  the 
Navier-Stokes  equation  in  spherical  co-ordinates  for 
spheres  of  Re=0.1  and  using  the  analytical  solution  of 
stokes  flow  for  spheres  of  Re  0.1  as  suggested  by  Le  Clair 
et  al  (1970).  The  hydrodynamic  interaction  between  an 
ice  crystal  and  a  water  drop  was  simulated  by  using  a 
superposition  model  of  their  combined  flow  fields.  The 
trajectories  of  water  droplets  approaching  ice  spheroids 
were  calculated  by  equations  of  motion  which  considered 
the  gravitational  force  and  the  drag  force  on  the  bodies 
under  consideration.  These  equations  were  integrated  by 
the  Hamming  method  and  initiated  by  the  Runge  Kutta  method. 
The  maximum  angle  that  water  droplets  would  collect  on 
the  edge  of  an  ice  crystal  was  also  investigated.  The 
flow  fields  around  an  oblate  spheroid  of  axis  ratio  0.05 
agreed  quite  well  v/ith  those  obtained  by  Bitter  et  al  (1973) 
but  the  trajectories  generated  in  the  present  study  were 
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significantly  different  from  those  of  fitter  et  al  (1974) 
in  that  droplets  approaching  the  center  of  the  spheroid 
were  not  forced  out  and  past  the  spheroid.  The  maximum 
angle  that  water  drops  accreted  on  the  edge  of  an  ice 
crystal  supported  the  theory  fo  graupel  growth  suggested 
by  List.  Hence,  a  thin  ohlate  spheroid  was  used  to  approxi¬ 
mate  a  hexagonal  ice  crystal  and  the  theory  of  Knight  and 
Knight  which  involved  corner  effects  on  accretion  patterns 
was  not  covered. 
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CHAPTER  I  INTRODUCTION  AND  LITERATURE  REVIEWS 


1 • 1  Objective  of  this  dissertation 


It  is  believed  that  hailsones  initiate  from 
smaller  pellets  called  ’graupel'.  According  to  Fletcher 
(1969)5  hail  represents  an  extreme  case  of  cloud  droplet 
accretion.  There  is  also  a  qualitative  difference  between 
a  hailstone  and  a  graupel,  A  hailstone  is  formed  by  a 
fast  accretion  process  during  which  the  heat  conducted 
away  from  the  stone  is  not  sufficient  to  freeze  all  the 
accreted  water.  Therefore,  part  of  the  water  remains 
unfrozen  and  the  stone  becomes  covered  w ith  a  liquid  skin 
which  freezes  relatively  slowly  to  a  dense,  more  or  less 
transparent  mass  with  air  bubbles  in  it.  On  the  other 
hand,  when  the  accretion  rate  is  slow,  each  collected 
cloud  droplet  freezes  almost  immediately  upon  impact.  The 
resulting  structure  is  porous  and  of  low  density  and  the 
particle  formed  is  identified  as  a  graupel  or  snow  pellet. 
Since  the  Nineteenth  Century,  there  have  been  many 
theories  involving  the  formation  of  graupel.  These  theories 
will  be  reviewed  in  the  later  section.  For  almost  two 
hundred  years  no  widely  accepted  formulation  of  the 
development  of  graupel  has  ever  been  provided.  The 
commonness  of  graupel  in  convective  storms  during  the  cool 
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days  of  Spring  suggests  that  they  are  one  of  the  fastest- 
growing  of  all  the  precipitation  forms  and  that  they 
must  be  common  within  summer  thundershower  clouds  although 

on  warm  summer  days  they  would  melt  into  raindrops  before 

% 

reaching  the  ground.  They  must  also  act  as  hailstone 
embryos  in  the  cold  regions  of  thunderclouds.  As  such, 
they  are  very  important  in  explaining  the  effectiveness 
of  rain-making  or  hail  suppression  mechanisms. 

There  is  no  strong  reason  to  suppose  that  there 
are  differences  in  origin  between  graupel  of  various  shapes. 
Dissection  of  small  hexagonal  graupel  always  clearly  reveals 
an  ice  crystal  within.  It  is  logical  therefore  to  believe 
that  graupel  start  as  a  rimed  ice  crystal  in  some  way  or 
another.  As  the  ice  crystal  continues  to  rime  it  could 
either  break  up  to  form  several  graupel  embryos  (Knight 
and  Knight,  1973)  or  it  could  continue  to  grow  into  a 
single  graupel  with  an  ice  crystal  at  its  center  (List, 

1958)  which  would  be  difficult  to  identify  in  a  laboratory 
dissection.  Since  the  shape  of  large  graupel  is  often 
conical  it  would  be  of  interest  to  know  the  side  of  the 
original  ice  crystal  on  which  the  riming  originally  begins, 
i,e,  fanning  out  into  the  wind  (Reynolds,  1876)  to  form 
the  blunt  end  of  the  cone  or  pointing  into  the  wind  (Aren- 
berg,  19Vl)  to  start  the  pointed  end  of  the  cone.  The 
main  interest  of  this  thesis  is  to  find  the  direction  of 
growth  of  a  graupel  pellet  during  its  early  stage  of 
growth  by  using  a  numerical  model  which  can  also  calculate 
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the  angle  of  the  point  of  graupel  which,  from  reports 
in  the  literature,  varys  from  60°  to  120°  (Arenberg,  1941; 
Flogel,  1877).  Wilkins  and  Auer  (1970)  reported  that  the 
onset  of  riming  is  mainly  on  plane  ice  crystals  between  300 
microns  and  800  microns  in  diameter.  These  atmospheric 
ice  crystals  typically  have  Reynolds  numbers  less  than  100 
and  greater  than  1  which  lie  between  the  values  described 
by  Stokes  and  potential  theories.  It  is,  therefore, 
necessary  to  obtain  solutions  for  the  viscous  flow  past  ice 
crystals  at  intermediate  Reynolds  numbers  in  order  to 
calculate  reasonable  flow  fields  and  hydrodynamic  interac¬ 
tions  between  ice  crystal  and  approaching  water  drops. 

1.2  Graupel;  its  definitions  and  seasonal  occurrence 

The  word  graupel  is  frequently  used  by  cloud 
physicists  and  its  definition  is  still  vague.  Evidently 
it  is  derived  from  the.  German  word  'Graupeln'  for 
hulled  barley  or  hulled  grain.  It  is  usually  translated 
in  German-English  dictionaries  as  ’sleet’,  an  ill-defined 
term  in  the  English  language  with  various  regional  interpre¬ 
tations,  or  'ice  pellet'  which  is  clearly  contrary  to  'snow 
pellets'  or  'soft  hail’  as  graupel  is  defined  by  Mason  (1971)* 
According  to  Mason,  the  International  Commission  for  Ice 
and  Snow  recommended,  in  19 56,  the  definition  of  graupel 
should  be:  "Snow  pellets  (soft  hail,  graupel).  These  are 
white,  opaque,  rounded  or  conical  pellets  of  diameter  up  to 
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about  6  mm.  They  are  composed  largely  of  small  cloud 
droplets  individually  frozen  together,  have  a  lou  density, 
and  are  readily  crushed.  They  may  break  up  on  striking  a 
hard  surface." 

According  to  Nakaya  (1951)  195*+)5  Schaefer 
(1951)  and  Magono  and  Lee  (1966),  graupel  are  the  extreme 
stage  of  riming  of  ice  particles  (usually  snow  crystals) 
before  further  accretion  leads  to  harder  growth,  making 
hail.  Some  hailstone  embryos  and  small  hailstones  are 
most  probably  soaked  or  partially  melted  and  then  refrozen 
graupel  (List,  1958). 

There  are  three  typical  forms  of  graupel: 
hexagonal,  conical  and  lumpy  or  granular  (Nakaya,  1951; 
Schaefer,  1951;  Magono  and  Lee,  1966).  Hexagonal  graupel 
are  obviously  heavily  rimed  snow  or  ice  crystals  (rimed 
plates  and  stellar  dendrites)  which  retain  their  general 
shape.  Conical  graupel  are  sometimes  pointed  and  sometimes 
not  and  they  occur  more  frequently  than  the  other  two  types 
in  convective  thunderstorms.  Lumpy  graupel  may  be  rimed 
snow  crystals  or  ice  crystals  or  parts  thereof  which  tumble 
along  in  the  process  of  riming. 

Very  often,  the  words  snow  crystal  and  ice 
crystal  are  used  interchangeably.  In  fact,  there  is  a 
difference  between  the  two  terms.  Ice  crystal  refers  to 
all  ice  particles  grown  in  an  ice-supersaturated  atmosphere 
by  diffusion  of  water  vapour  to  their  growing  surfaces,  and 
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may  exist  as  individual  units  of  simple  geometric  shape, 
for  example  hexagonal  prismatic  columns  (needles  and  prisms) 
and  hexagonal  plates,  or,  under  suitable  conditions,  may 
grow  into  complex  richly-branched  forms.  Snow  crystal 
includes  all  solid  precipitation  in  nature.  They  can  grow 
by  any  of  the  three  possible  process:  (1)  growth  by 
diffusion  in  the  presence  of  supercooled  drops,  (2)  growth 
by  collision  with  other  solid  hydrometeors  in  the  atmosphere, 
or  (3)  growth  by  coalescence  of  supercooled  water  droplets. 
Thus,  the  term  snow  crystal  can  be  implied  for  ice  crystal, 
but  not  vice  versa.  The  above  definitions  of  these  two 
terms  will  be  used  in  this  study. 

Graupel  fall  to  the  ground  during  three  rather 
distinct  seasonal  weather  phenomena.  These  are  (1),  heavy 
winter  convective  snow  storms  (2),  convective  storms  in 
the  early  Spring;  and  (3),  as  hailstone  embryos  in  Summer. 

It  is  rare  to  find  hexagonal  graupel  in  vigorous  convective 
storms  or  occuring  as  hailstone  embryo,  (Knight  and  Knight, 
1973)  and  conical  graupel  predominates  in  all  three  of  the 
above  storms. 

1.3  Review  of  the  theories  of  conical  graupel  formation. 

The  origin  of  the  three  seasonal  occurrences  of 
graupel  mentioned  above  is  probably  quite  similar.  Every 
investigator  agrees  and  the  observational  evidence  shows 
that  graupel  are  a  result  of  riming;  i.e.  accretion  of 
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supercooled  water  droplets  in  conditions  such  that  each 
drop  retains  its  shape  when  it  freezes  on  impact.  The 
droplet  becomes  a  little  ball  of  ice  adhering  to  the  object 
being  rimed. 

While  the  initiation  of  snow  pellets,  or  graupel, 
is  generally  understood  to  be  due  to  the  riming  of  settling 
snow  crystals,  there  are  certain  features  in  their 
development,  particularly  the  conical  form,  which  have  not 
achieved  a  widely  acceptable  description. 

Arenberg  (19Vl)  proposed  that  conical  graupel 
begin  as  large  snow  crystals  that  become  heavily  rimed, 
on  being  carried  upwards  to  the  top  of  cumulus  cloud.  If 
the  flakes  are  flat  hexagonal  plate  of  various  degrees  of 
intricacy  they  would  fall  with  their  flat  side  facing  down¬ 
ward  vertically  as  rime  continues  to  collect  on  their  under 
side.  First  they  would  rime  along  the  edges  and  then  over 
the  entire  surfaces.  The  deposit  gradually  builds  up  until 
the  point  of  the  cone  is  formed.  Arenberg  admits  that  this 
is  not  the  most  stable  fall  attitude  for  a  cone,  i.e.  pointed 
side  down,  but  he  points  out  that  it  is  the  second  most 
stable  attitude.  It  could  be  assumed  that  during  the  early 
growth  stages  when  there  is  little  difference  in  stability 
between  the  two  stable  fall  modes  that  the  pointed  downward 
mode  could  be  maintained  for  long  periods  if  no  disturbance 
occurs.  His  proposition  explains  also  the  observation  that 
graupel  are  least  dense  at  the  top  and  harder  at  the  base. 
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Holroyd  (196*+)  suggests  that  conical  graupel 
start  as  rimed  clusters  of  two  to  six  rimed  needle  crystals, 
as  deduced  from  observations.  He  proposed  that  needles 
falling  in  their  usual  horizontal  orientation  would  tilt 
and  descend  vertically  as  rime  build-up  makes  one  end 
heavier  than  the  other.  The  increasing  deposits  of  rime  on 
the  downward  end  keeps  the  needles  in  that  orientation.  If 
two  rimed  needles  coalesce  they  would  attach  at  the  bases 
where  their  dimensions  are  greatest  and  the  air  flowing 
past  would  then  force  the  trailing  tips  together.  Further 
riming  or  more  coalescence  with  other  rimed  needles  would 
develop  a  conical  shape,  A  conical  graupel  would  have  been 
formed. 

While  most  investigators  believe  that  conical 
graupel  start  as  ice  crystals,  studies  by  Weickmann  (1953) 
and  Magono  (1953)  suggested  that  conical  graupel  can  start 
directly  from  frozen  droplets  or  tiny,  irregular  ice 
particles.  Their  proposition  was  based  mainly  on  the  fact 
that  snow  crystals  are  seldom  recognized  within  conical  hail 
embryos. 

Contrary  to  Arenberg*s  pointing  downward  theory, 
most  writers  since  the  very  early  investigation  of  graupel 
by  Reynolds  ( 1 876 )  expect  riming  to  lead  to  shapes  that  grow 
broader  on  the  bottom  side,  i.e.  toward  the  approaching  air 
and  cloud  droplets.  In  such  a  growth  mode,  the  original 
snow  crystals  must  then  be  near  the  tips  of  the  graupel, 


\ 


J 


8 


not  at  the  blunt  ends.  List  (1958)  shows  photographs  to 
illustrate  that  ice  crystals  can  be  found  near  the  pointed 
ends  and  he  (I960)  also  points  out  that  about  h-0  or  50 
percent  of  the  graupel  he  examined  in  Switzerland  had 
originating  snow  crystals  as  centers.  If  this  is  the 
actual  way  of  most  graupel  growth,  for  the  pointed  ends  to 
be  truly  pointed  some  capture  must  occur  on  the  down  wind 
side  of  the  original  snow  crystal,  i.e.  in  the  wake. 

Most  recently,  Knight  and  Knight  (1973) 
hypothesize  that  conical  graupel  should  form  as  follows. 

A  planar  or  dendritic  ice  crystal,  falling  in  its  most 
stable  orientation  with  its  axis  vertical,  tends  to  rime 
first  at  the  edges  of  the  windward  side,  in  particular  it 
rimes  more  at  the  six  corners.  The  riming  starts  at  points 
of  the  crystal  and  fans  out  into  the  wind.  Further  riming 
on  the  crystal  which  is  still  falling  in  the  same  attitude 
increases  the  size  of  the  conical  shape  deposits  attached 
by  their  tips  to  each  of  the  six  corners  of  the  snow  crystal. 
If  turbulence  exists  or  the  snow  crystal  changes  its 
orientation,  these  conical  deposits  would  break  off  due  to 
air  stress.  Graupel  could  then  originate  as  single  cones, 
which  have  broken  off  from  the  corners  of  ice  crystals,  or 
as  groups  of  cones  that  have  grown  together  and  coalesced 
at  their  broader  ends.  This  would  explain  the  difficulty 
in  finding  ice  crystals  within  graupel  particles. 

The  fact  that  the  collection  efficiency  is  the 
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greatest  at  the  edges  of  the  downward  facing  side  of  a 
falling  snow  crystal  of  size  1  50  to  400  microns  has  been 
justified  by  Pitter  and  Pruppacher  (197*+)  through  numerical 
investigation.  Sasyo  (1971) »  after  careful  study  of  the 
formation  of  precipitation  by  the  aggregation  of  snow 
particles  and  the  accretion  of  cloud  droplets  on  snowfakes, 
also  provided  evidence  that  cloud  droplets  coalesce  in 
greater  amounts  at  the  corners  than  over  the  other  parts  of 
snow  crystal.  Using  supported  models  of  hexagonal  plates 
of  diameter  0.5  cm  and  1.0  cm  placed  in  saturated  air  moving 
at  various  relative  velocities,  he  found  that  the  frontal 
surface  captures  water  droplets  of  various  sizes  while  the 
rear  surface  captures  only  water  droplets  smaller  than  10 
microns  in  diameter.  Furthermore,  he  discovered  that  there 
is  a  tendency  to  increase  the  number  of  droplets  captured 
by  the  frontal  surface  with  increasing  air  velocity  at  the 
same  time  as  droplets  captured  by  the  rear  surface  increased. 
It  seems  that  Knight  and  Knight’s  theory  of  graupel  growth 
is  supported  by  these  studies  of  Sasyo. 


CHAPTER  II  VISCOUS  FLOW  PAST  ICE  CRYSTALS  AND  CLOUD  DROPS 


2.1  Introduction 


The  first  step  in  this  numerical  study  is  to 
calculate  flow  fields  around  ice  crystals  and  cloud 
droplets.  The  second  step  is  to  calculate  the  motion  of  a 
cloud  droplet  as  it  approaches  an  ice  crystal.  There  are 
numerous  reports  of  experimental  and  theoretical  work  in 
the  literature  on  the  characteristics  of  the  flow  of  a 
viscous  fluid  past  a  single  sphere  and  on  the  forces  between 
two  hydrodynamically  interacting  spheres.  (Le  Clair  et  al, 
1970,  1972;  Hamielec  et  al,  1967;  Shafrir  et  al,  1970). 

In  contrast,  no  work  has  been  done  on  hexagonal  discs  and 
relatively  little  on  planar  objects.  Jayaweera  and  Cottis 
(1969))  and  List  and  Schemenauer  (1971)  experimentally 
showed  that  the  hydrodynamic  behaviour  of  a  simple  hexagonal 
plate  can  be  approximated  with  sufficient  accuracy  by  that 
of  a  circular  disk.  Pitter  et  al  (1973)  showed  that  the 
hydrodynamic  behaviour  of  a  thin  oblate  spheroid  is  essential¬ 
ly  the  same  as  that  of  a  circular  disk  of  the  same  thickness 
to  width  ratio  (aspect  ratio).  The  flow  field  past  a  disk 
can  most  easily  be  calculated  if  its  shape  is  approximated 
as  a  thin  oblate  shperoid  of  variable  aspect  ratio.  The 
flow  field  past  a  spherical  cloud  droplet  would  be 
equivalent  to  an  aspect  ratio  equalling  'one'.  The 
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interaction  of  the  calculated  flow  fields  around  the  plate 
and  the  sphere  should  allow  one  to  calculate  the  trajectory 
of  a  cloud  droplet  as  it  is  approaches  and  collides  with 
an  ice  crystal. 

2.2  The  equations  representing  flow  around  a  spheroid 
or  sphere 

The  motion  of  a  system  can  be  best  described 
by  Newton’s  second  law  as: 

(mass  of  system) (acceleration  of  the  systera)=net  force  of 
the  system.  For  a  fluid  parcel,  the  net  force  can  be 
classified  as  the  body  (gravity,  etc.)  force  and  the 
surface  force.  Newton’s  second  law  can  then  be  expressed 
in  an  equational  form  as, 

DV  net  body  force  net  surface  force 

density  —  =  -  +  -  , 

Dt  volume  volume 

where  the  time  derivative  of  the  velocity  is  the  acceleration 
.DV 

(—  )  of  the  fluid  parcel.  For  a  Stokesian  fluid,  such  as 
air,  water  and  oil,  the  net  surface  force  on  a  fluid  parcel 
per  unit  volume  is  formed  by  two  forces;  the  pressure 
gradient  force  and  the  viscous  friction.  If  the  net 
body  force  per  unit  volume  is  denoted  by  F,  and  the  fluid 
is  incompressible,  then  the  equation  of  motion  and  continuity 
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equation  can  be  written  as, 


DV*  9  _ 

p- -  =  F  -  VP*  +  M*V^V*  (2.1) 

Dt  ' 


V*V*  =  0  (2.2) 


where  P*  and  V*  are  the  dimensional  pressure  and  velocity, 

p  and  yx*  are  the  dimensional  density  and  dynamic  viscosity 

of  the  fluid.  The  differential  operators:  gradient  V, 

2 

divergence  V*  and  Laplacian  V  ,  in  generalized  curvilinear 
orthogonal  co-ordinates  are: 


V  -  +  S?£-  -t  5_3^- 

h,  h2  Bli  , 


1 


VV  = 


h1h2h3  . 


k  *1 


3q. 


(V3h1h2) 


/ 


(2.3) 

(V2h3h-} )  + 

(2.1) 


and 


1 


h1h2h3 


3 _  (  h3h1 

£q2  h2 


^  h2h3  J) _ 

h1  c)  q1 


+ 


)  +1_  (Hi 

3  q2  3q3  h3 


(2.5) 

Q-j  »  ^2  ?  ^3  ’  are  curv^-^^near  co-ordinates; 

3^  ,  ai  ,  a3  are  unit  direction  vectors  and  h^ ,  h2,  h^ 
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are  the  metric  coefficients.  The  values  of  the  metric 
coefficients  used  in  the  co-ordinate  system  of  this 
study  are  given  in  section  2.*+.1.  If  the  fluid  is  steady 
state  equation  (2.1)  becomes: 

V*.VV*-yuL*V2  V=  -  VP  +  ?  (2.6) 


The  calculation  of  the  flow  fields  in  this 

study  is  restricted  to  axial  symmetry,  i.e.  there  are  no 

variations  with  the  azimuthal  angle.  A  three  dimensional 

approach  would  be  much  more  complicated.  Once  axial 

symmetry  is  imposed,  i.e.  =  0  and  =  constant,  one 

3  1 1 1  * 

can  define  a  Stokes’  stream  function  t  which  will  satisfy 
the  continuity  equation  (2.2)  automatically,  and  use  it  as 
one  of  the  dependent  variables  to  eliminate  the  pressure 
gradient  term.  If  the  primitive  variables,  velocity  and 
pressure,  are  used;  special  care  must  be  taken  to  satisfy 
the  continuity  equation  as  accurately  as  possible.  The 
stream  function  t^*can  be  defined  in  terms  of  velocity 
in  the  two  dimensional  curvilinear  co-ordinate  system 
with  q-|  as  the  radial  and  as  th0’  angular  coordinates: 


V* 


V* 


-i  ay' 


h2h3  3  q2 


1 


h1h3  aq, 


(2.7) 
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The  other  dependent  variable  is  the  vorticity  fl*  =  V  x  V*, 
Because  of  the  symmetry,  the  only  velocity  component  that 
survives  is  the  one  in  the  azimuthal  direction.  This 
component  is  designated  by  ^  • 

If  a  curl  operator  is  applied  to  equation  (2,6) 
the  pressure  gradient  term  becomes  zero,  i.e,  V  x  VP  =  0 
and  the  operation  of  the  curl  on  the  velocity  vector  in  a 
axisymmetric  curvilinear  coordinates  gives: 


V  x  V*  = 


h,  h 


SluV,,  3h,V„* 


z  ^ 


3  4' 


aq; 


A  'yS 


(2.8) 


then  equation  (2,6)  becomes: 


where 


i  aO*,  e*VA*) 
a  c  q, ,  qz ) 

d  (  A  .  a  t  h, 


(2.9) 


h.hj 


*s 


|V  h. h3  dq,  '  ^2^3  ^  ^2 


■)! 

j  • 


The  operator  5  (  E*2-  %/ h 3  )  / d  (q^i)  is  a 

Jacobian  and  it  is  denoted  by  the  determinant  as  follows: 


a(^*,E*V/h3) 


9(E*V7hb 
3q,  aq2 

9(e*Wh?)  a^* 

aq,  aqz 


(2.10) 
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The  other  operator  E*  is  given  by: 


E* 


2  _  hs  d  /  hz  _ \  c)  /_Jh_  _c) _ ^ 

H4h2 L  d<3,\  h,h3dq/ 


a 


(2.11) 


The  symbol  V*  in  equation  (2.9)  is  the  kinematic  viscosity 
and  it  is  defined  as  the  dynamic  viscosity  yU*  divided  by 
the  density  p.  Equation  (2.9)  is  of  the  fourth  order  in 
Vp*  and  it  is  not  linear  because  of  the  presence  of  the 
inertial  term  on  the  right  hand  side.  In  general,  the 
equation  is  called  the  Navier-Stokes  equation  in  two 
dimensions. 


2.3  Analytic  Solutions  to  Flow  Around  Spheres  and 
Spheroid  s 


The  oldest  formulation  of  flow  past  an  obstacle 
is  called  Stokes  flow  (18^1)  or  creeping  flow.  The 
inertial  effects  are  assumed  to  be  negligible  in  comparison 
with  those  of  viscosity.  The  Navier-Stokes  equation 
reduces  to: 

E*2E*2lp*  =  0  (2.12) 

Such  a  simplification  of  the  equation  of  motion  is  very 
substantial,  since  the  equation  of  motion  thus  becomes 
linear. 
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The  general  solution  of  the  creeping  flow 


equation  (2.12)  for  spherical  (r,  8)  and  spheroidal 
coordinates  (  )  was  given  by  Sampson  (1891).  For 

spheres , 


+  Dnr-«+3)lnU) 


/  n+2 


(2.13) 


■where  An,  B  ,  Cn,  and  the  corresponding  primed  terms 
are  constantsand  Y  is  defined  as  Y  =  cos  0.  In(  Y  )  and 
Hn(  Y)  are  G-egenbauer  functions  of  order  n  and  degree  -i 
of  the  first  and  second  type  respectively.  The  constants 
are  evaluated  from  the  applicable  boundary  conditions. 

For  the  no-slip  condition  at  the  sphere  and  a  uniform 
stream  function  at  infinity,  equation  (2,13)  becomes: 


y  *  =  i  U*  a2  ( +  ~  +  )  sin2  0  ,  (2.14) 


■where  a  is  the  radius  of  the  sphere  and  U*  is  the  speed 
of  the  undisturbed  stream.  The  drag  force  experienced  by 


the  sphere  can  be  given  as: 


F  =  4  IT yU.*D2  , 


(2.1?) 


. 
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T>2  is  the  constant  from  equation  (2.13)°  It  is  equal  to 

—  a  U*.  The  dimensional  surface  pressure  is  given  by: 

2 


P*=  p  u*2  (  1  +  ~  cos  0  ).  (2.16) 

Re 


For  Stokes  flow  in  the  spheroidal  coordinates  (  )5 

\y*  was  given  by  Sampson  (1891)  and.  discussed  in  Masilyah 
(1973)  as: 

IV*  =  C0+  D.T+  +  C2I»(t)  +  D,Hxtt) 

-t- B4. H4 (x)^  +  [b5  +  D5h3(t;) -t  B5  H. 


+  B4H2(t)  +  D4H4(t)+  B,'H 


6'  % 


(2.1 7) 


where  X=cosht  ;  T=coS1j  InW  are  the  Gegenbauer 
functions  of  the  first  type  with  corresponding  orders; 

Hn(  )  are  the  Gegenbaver  functions  of  the  second  type 
with  corresponding  orders;  and  all  the  C's,  Brs  and  Ds  with 
or  without  primes  are  constant  coefficients  dependent  on 
boundary  conditions.  For  an  oblate  spheroid  one  obtains: 


q* *  =  +  t  (2.18) 

where  i  =  J -1 ,  the  imaginary  unit,  and  /\=  sinh ^ •  The 
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spheroidal  coordinates  are  given  as  in  figure  (2«1).  For 
the  no-slip  condition  at  the  spheroid  surface  and  a 
uniform  stream  condition  at  infinity,  the  stokes  stream 
function  given  by  equation  (2.18)  becomes: 


4>*  =  iu*c2  ( /\2+  0(  i  -'f 2) 

l  [/\a/(/\2a+0]-[(^a-0/(/\a+l)]cot"Az(J-i9) 

V  / 

where  subscript  OL  denotes  the  surface  and  c  the  focal 
length  of  the  spheroid. 

The  drag  force  experienced  by  the  spheroid  is 
then  given  as; 


F  =  6  TTyU-*  U*QeK  t  (2.20) 

where  CLeis  the  major  axis  radius  of  the  spheroid,  and  for 
an  oblate  spheroid, 


_ 1 _ _ 

I  £  /\a"  (Aa  “I ) c°t  1  ^ a 3  • 


(2.21) 


The  dimensionless  surface  pressure  distribution  for  Stokes 
flow  past  an  oblate  spheroid  is  given  by 


P  = 


I  -f 


6  K 


Xq  t  1  \ 

\a  *  Cosl^J 


Re 


(2.22) 
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In  the  case  where  the  inertial  effects  of  the 
flow  are  small  but  not  negligible,  Stokes  solution  can  no 
longer  be  considered  valid  at  points  far  away  from  the 
surface  unles  Re  0.  At  points  near  the  body  the  inertial 
forces  tend  to  vanish.  This  means  that  Stokes  approximation 
for  the  flow  field  is  not  uniformly  valid  throughout  as 
it  breaks  down  at  a  large  distance  from  the  surface. 

The  inertial  terms  are  thus,  to  some  extent,  taken  into 
account,  Goldstein  (1929)  developed  Oseen's  (1915) 
perturbation  method  of  approximating  the  drag  force 
experienced  by  a  sphere  in  non-Stokes  flow,  giving  the 
drag  force  as  a  series  in  Reynold  number: 


F  = 


61TyU*aU* 


1  + 


3_ 

16 


1 9  Re 

1280 


4* 


1 1  -Re! 

204-80 


3oil9Re4 

34406400 


122519  Re5 

56°74240 


#  •  •  ♦  * 


1 


(2,23) 


The  first  term  of  the  above  equation  is  the  value  given  by 
Stokes  (equation  2, 15),  while  the  second  term  is  that  given 
by  Oseen,  Unfortunately,  the  equation  (2,23)  is  valid  only 
for  Reynold  number  up  to  and  including  2, 

Aoi  (1 955) ,  maintaining  conformity  with  the 
corresponding  creeping  flow  solution,  derived  the 
expression  for  the  total  drag  force  on  an  oblate  spheroid 
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as: 


16 


(2.24) 


where  K  and  Cle  are  the  same  as  those  in  equation  (2.18), 
while  the  Re  is  based  on  the  diameter  of  the  major  axis 
of  the  spheroid.  Aoi  (19 55)  also  developed  the  ratio  of 
the  form  (pressure)  drag  to  the  skin  (friction)  drag  on 
an  oblate  spheroid  which  is  given  by: 


F P  (Xa  +  lK  I  -Xa cot'' Aa ) 


I - (\a  +  I ) ( 1  -  AqCoI  'Aa) 


(2.25) 


Thus,  according  to  the  above  expression,  though  the  form 
drag  and  skin  drag  of  the  oblate  spheroid  depend  on  its 
size,  the  ratio  of  the  drag  forces  is  independent  of  the 
Reynolds  number. 


Proudman  and  Pearson  (1957)*  instead  of 


obtaining  perturbation  fields,  seek  separate  asymptotic 
solutions  which  are  locally  valid  in  separate  regions 
near  to,  and  far  from,  the  body.  These  inner  and  outer 
solutions  are  each  determined  by  asymptotically  matching 
them  in  their  common  domain  of  validity.  The  expression 
for  the  drag  force  experienced  by  a  sphere  placed  in  a 
uniform  stream  is  given  by: 


(2.26) 


- 
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Breach  (1961)  extended  the  technique  of  Proudraan  and 
Pearson  to  apply  to  ellipsoids  of  revolution,  both  prolate 
and  oblate.  The  expression  for  the  drag  force  on  an  oblate 
spheroid  can  then  be  rewritten  in  the  form: 

F  =  6TtAA*U*aeK(  i  +  lRe  +  3.  K2Re  In  Re2)  (2,27) 

7  l  16  160  Z  )  y 

with  Qe  and  K  identical  to  those  in  equation  (2,20),  Once 
again,  the  first  term  in  the  above  equation  is  for  creeping 
flow,  the  second  being  Oseen's  extension  while  the  third 
is  new.  For  a  sphere,  K=1  and  the  equation  (2,27)  becomes 
identical  to  that  of  Proudman  and  Pearson  (equation  2,26), 

2.4  Numerical  Solutions  of  Flow  Around  Spheres  and 
Spheroids. 

The  solution  of  the  Navier-Stokes  equation  via 
Stokes,  Oseen  and  Proudman  and  Pearson  fail  to  describe  the 
flow  for  Reynolds  number  greater  than  2,  Approximate 
solutions  using  the  Galerkin  method  (Snyder,  Spriggs  and 
Stewart,  1964)  have  been  tried  by  numerous  investigators. 

This  method  first  assumes  trial  functions,  and  an  approximate 
solution  is  then  obtained  by  determining  the  trial  function 
parameters  such  as  to  satisfy  the  Navier-Stokes  equation 
and  the  boundary  conditions  as  closely  as  possible.  The 
optimum  trial  functions  may  be  determined  by  variational 
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calculus  or  by  error  distribution  methods.  The  success  of 
this  method  depends  strongly  on  the  choice  of  the  trial 
polynomial,  Hamielec  and  Johnson  (1962),  Hamielec  et  al 
(1963)  used  thus  technique  to  evaluate  the  viscous  flow 
around  fluid  spheres  of  Reynolds  number  from  10  to  5000, 

Jenson  (1959) 5  in  his  work  on  spheres,  solved 
the  Navier-Stokes  equation  by  splitting  it  into  two 
simultaneous  second  order  equations  representing  the  stream 
function  and  the  vorticity,  Jenson  approximated  the 
vorticity  near  the  surface  by  a  third  order  polynomial  and 
he  used  an  exponential  transformation  for  radial  distance 
from  the  surface,  Le  Clair  and  Hamielec  (1970)  and 
Hamielec  et  al  (1967)  extended  Jenson's  technique  with  the 
aid  of  digital  computers  to  intermediate  and  higher 
Reynolds  number, 

Rimon  and  Lugt  (1967),  using  a  time-dependent 
numerical  technique,  calculated  the  viscous  flow  around 
an  oblate  spheroid.  Unfortunately,  their  results  were 
limited  to  Re  =  10  and  Re  =  100,  Masliyah  and  Epstein 
(1971)  using  an  adaptation  of  the  relaxation  technique 
of  Jenson,  obtained  numerical  solutions  of  the  Navier- 
Stokes  equation  for  axisymmetric  flow  past  oblate  spheroids 
at  particular  Reynolds  numbers  up  to  100  and  of  aspect 
ratios  0o2  to  0,999«  Most  recently,  Pitter  et  al  (1973) 
followed  the  numerical  method  of  Woo  (1971),  adapted  to 
oblate  spheroids  and  obtained  quite  satisfactory  solutions 
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of  the  steady-state  Navier-Stokes  equation  of  viscous  flow 
past  an  oblate  spheroidal  obstacle  of  Reynolds  numbers  0.1 
to  100  and  of  axis  ratio  0.05  to  0.2.  Their  numerical 
results  for  oblate  spheroids  of  aspect  ratio  0.2  agree 
well  with  those  of  Masliyah  and  Epstein  (1971)* 

2.5  Finite  Difference  Equation  Formulation  of  Plow  Around 

Spheres,  and  Spheroid, 

The  choice  of  spherical  (for  water  drops)  and 
spheroidal  (for  ice  crystals)  coordinates  in  this  study 
facilitates  the  formulation  of  the  boundary  conditions.  In 
addition,  if  these  coordinates  are  given  exponential 
properties,  a  fine  lattice  near  the  surface  of  the  particle 
and  a  coarse  lattice  far  away  from  the  body  can  be  formulated. 
As  the  influence  of  the  particle  on  the  flow  is  mainly 
manifested  near  the  particle's  surface  and  as  this  influence 
decreases  with  distance  from  the  surface,  the  lattice  thus 
produced  serve  the  purpose  of  giving  a  desired  lattice 
gradation. 

2.5*1  Spheroids 

The  Navier-Stokes  equation  for  steady-state 
axisymmetric  or  two-dimensional  flow,  expressed  in  orthogonal 
curvilinear  coordinates  in  terms  of  the  dimensional  stream 
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function  Lp  *  is: 


V  *  E*2  E*2  \\)  * 


1 

ti  -j  H^Ti^ 


3CV*,  E*^*/h3*) 

acq,,qa) 


(2.9) 


and 


2  _  h3 


E*  = 


h,h2 


B  f  h2  3 _ ^  L~  h _ (  K  A. 


■) 


dQph.ha,  5^/  ^2.1^5  5  (\i  J  . 


(2.11) 


In  spheroidal  coordinates, 


q1  =  ? 

*2  =  4 


h^  =  h2  =  c  (Sinh  ^  +  cos^  )•§■ 


h 


^  =  c  cosh  ^  sin  y|^ 


where  ^  is  the  radial  and  ^  is  the  angular  coordinates  as 
in  figure  (2.1)  and  C  is  the  focal  length  of  the  system, 
i.e.  the  distance  between  the  two  foci  of  the  oblate 
spheroid,  then,  equation  (2.9)  and  equation  (2.11)  become: 


y  *  E*2  E*2  Vj/  * 


and 


(  cosh  E; 

CiCs(nh2^  +  cosi^) 

e 


(2.28) 


E*2  = 


CJ(si«h2§  +co^ii)l.  3P 


["— - tanli£".iL-  +  3 _ Cot  ii  3 
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Define  the  dimensionless  stream  function 
=  4<yu*  a1  ;  the  dimensionless  Reynolds  number  Re= 
2aU^V*  ;  and  set  E^  =  a^E*^  where  ’a'  is  the  semi¬ 
major  axis,  a  =  C  cosh  ^a.  ^a  is  the  value  of  5;  at  the 
surface.  Equation  (2*28)  and  Equation  (2*29)  become: 


Sech  £aE  7(  E 


Re  I 

Cosh£;  SU\>| 

z 

+  Cos**} 

f  sty  5  / 

E2ty 

\  d  (  E2tp  \ 

^  d^l 

cosh^sln  ^ 

j  d^cosh^sin  ^ 

(2. 30) 


Cosll2.?  q 


SvnH 


.tanh  6  A 


2  * 
3 _ cot  v)  jL_ 


(2.31 ) 


Putting  E21|<  =  ^  cosh  sin  sech^^  ,  where  ^  is  the 

non-dimensional  vorticity  and  and 

2 

E*  V.  is  the  value  of  ^  at  the  surface,  then 
equation  (2*30)  can  be  split  into  two  equations, 


sech  £;a  EZ(  ^ cosh ^ s\n  W )  =  _Re  cosh  sm  ^ 

2.  ls(nk2^-t  co^r^j 

S 


5  1 

( cosh^s<*\iq 

(2*32) 
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Sinh2^  +  cos2^ls 

r=  ^  cosh  ^  scvq  ^  sech5 
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Since  ^  can  become  very  large  at  the  surface  of  the  spheroid 
and  hard  to  be  handled  numerically,  it  is  more  convenient 
to  transform  the  equation  (2*32)  and  equation  (2,33)  to 


E2G  = 


Re  cosh^ 
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r 


coskFsin.  Y[ 


+  Co 

^  2>F  3^  aF 


3^  3^ 


=  G  seek" 


a 


(2.31+) 


(2.35) 


by  defining  the  non-dimensional  quantities: 

(2.36) 

(2.37) 

Using  the  finite  differancing  scheme: 

_  d  $  j  __  jj-H  ~  St-' 
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where  A  and  B  are  the  radial  and  angular  spacing  of  a 
spheroidal  mesh  as  shown  in  figure  (2.2),  equations  (2.3^-) 
and  (2*35)  can  be  converted  into  finite  difference  form  as 
given  by  Masliyah  (1973)** 


2  2 

G(l , J)  (  — - -  +  — ^r)  =  G(I,J+1) 
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Equations  (2.38)  and  (2.39)  can  then  be  solved  simultaneously 
with  appropriate  boundary  conditions  of  the  dimensionless 
stream  function  and  non-dimensional  vorticity  by  computer 
to  obtain  corresponding  values  of  S'  and  ^  at  each  grid  point. 

2.5*2  Spheres 


The  non-dimensional  Navier-Stokes  equation 


for  steady-state  axisymmetric  flow  in  spherical  co-ordinates 
(r,9)  can  be  expressed  as: 


where 


(2.41) 


Putting: 


(2.42) 


where  iS  as  before  is  the  dimensionless  vorticity,  then 
equation  (2.40)  can  be  reduced  to: 


(2.43) 


and  equation  (2.42)  becomes: 


(2.44) 
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For  the  sake  of  having  exponential  properties  which  would 
give  a  fine  lattice  near  the  surface,  the  radial  coordinate 
r  is  transformed  into  r  =  ez.  The  equation  (2*43)  and 
(2*44)  are  now  given  as, 


R, 


3'4)  9 
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^  \  W  3  I  <  \~l 
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As  for  the  same  reasons  mentioned  in  the  case  of  spheroid, 
the  transformations: 

F  =  ^  /  ez  sin©  (2*47) 

G  =  C,  ez  sin©  (2.48) 


are  taken.  The  Navier-Stokes  equation  in  spherical 
co-ordinates  is  finally  expressed  as: 
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Taking  the  finite  differencing  scheme  as  before  with  C  and  D 
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as  the  radial  and  angular  spacing  respectively  of  a  spherical 


mesh  shown  in  fig*  (2*3)5  the  equations  (2. *+9)  and  (2*50), 
expressed  in  finite  differencing  form,  can  be  written  as: 


(  F  ( l+1,T)-F(i-un-  ( Wl+ljW.M  J)VF(  I, MV  F(I, 
^  2D  '  2D  A  2C 


(2.52) 


The  equations  (2*51)  and  (2*52)  can  be  solved  simultaneously 
by  computer  for  appropriate  boundary  conditions* 


2*5*3  Formulation  of  Boundary  Conditions  for  Spheroids. 

The  numerical  analysis  is  confirmed  between  an 
inner  surface  represented  by  the  spheroid  itself  and  by 
an  outer  envelope  which,  like  the  spheroidal  surface, 
coincides  with  one  of  the  spheroidal  co-ordinate  grid  lines. 

As  the  two  flow  equations  are  of  second  order, 
four  boundary  conditions  must  be  satisfied.  The  boundary 
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conditions  for  the  stream  function  ip  are: 


for  n  =  0,  ^  =0 

for  Y\jz  71,  <\)  =0 

for  ?  =  »  ^=0 


defines  the  axis  of  symmetry, 


at  the  surface  of  the  spheroid* 


At  the  outer  boundary,  the  flow  is  assumed  to  be  a 

p  2  f* 

streaming  parallel  flow,  giving  =  ■§•  sin  ^  cosh 
sech^J^'^  *  The  boundary  conditions  for  the  vorticity  ^  are 


for  =0 

for  =  Tf 


5  =0 
$  =0 


defines  the  axis  of  symmetry, 


for  ^  , 

f  0r  £  =  %co  > 


=E^  '-V  cosh^?  /sinfl  —  at  the 

t  surface  of 
the  spheroid 

=0  -  in  the  undisturbed  flow. 


The  boundary  condition  for  ^  at  ^  ,  which  originates 

from  equation  (2*35) 5  can  be  explained  as  follows:  At  the 
surface,  the  no-slip  conditions  noted  above  require  that 

ip  -  =  yiy  =  n 

Equation  (2.35)  therefore  reduces  to 
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J=1  is  the  value  of  J  at  the  surface  of  the  spheroid. 
Using  the  values  of  ^  at  J=2  and  J=3  and  expanding 
equation  (2,53)  as  a  Taylor’s  series  to  the  third  order, 
gives 


G(I,1) 


cosh2 8  Vd, 2)^(1, 5)J 

2  [* S^nh2  +  cos2  Kj(L)J 

^  9 


which  is  equivalent  to 


(2.5V) 


cosh??a[8<4JU,2)-Hni,3f 
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2.5«V  Formulation  of  Boundary  Conditions  for  Spheres. 


The  no-slip  boundary  conditions  for  a  sphere 
are  about  the  same  as  those  for  spheroid,  namely  for 
they  are: 

for  0  =  0  ,  \\)  =  0 

for  0  =  Tf  ,  0 

r  =  1  ,  4>=  0  _ 


Far  from  the  sphere,  there  is  undisturbed  parallel  flow: 
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The  boundary  conditions  for  are: 


for  9  = 

:0 

5 

5  = 

0 

0  = 
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E^  ^ /sinQ 
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■V 
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• 

The  boundary 

conditions 

for 

^  at  r=1 ,  can  be  as  follows* 

the  surface,  the  no-slip  conditions  require: 

Ul  -IV  -  -  0 

Y  36  ar  902 

Then  equation  (2,IP+)  reduces  to: 

_|ll  =  ^rsme  (2.5 

Expanding  as  a  Taylor’s  series  to  the  third  term,  equation 
(2*56)  takes  the  following  finite  difference  form: 

Cm  i)=  8^d,2)-im,5) 
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CHAPTER  III  HYDRODYNAMIC  INTERACTION  OF  AN  ICE  CRYSTAL 


AND  A  CLOUD  DROPLET 

3*1  Introduction 

The  basic  shape  of  ice  crystals  grown  by 
diffusion  is  that  of  a  hexagonal  prism*  Depending  on  the 
ambient  temperature,  growing  ice  crystals  possess  either 
prominent  c-axis  growth  (hexagonal  columns)  or  a-axis 
growth  (thin  .hexagonal  plates).  In  atmospheric  clouds 
the  warmest  temperature  at  which  glaciation  begins  is 
generally  between  -10°to  -18°C  giving  ice  crystals  with 
predominantly  plate  type  growth0  That  temperature  range 
corresponds  to  growth  mainly  in  the  form  of  thin  hexagonal 
plates  with  or  without  dendritic  extensions.  The  dendritic 
(i.e,  intricate)  crystals  are  caused  by  fast  growth  due  to 
high  supersaturations.  These  facts  have  incited  the  study 
of  the  collision  behaviour  of  plate-like  ice  crystals 
with  water  drops  as  a  study  of  the  riming  behaviour  of 
ice  crystals,  Pitter  and  Pruppacher  (197*0  computed  the 
collision  efficiencies  of  water  spheres  with  oblate  ice 
spheroids  of  axis  ratio  0,05  and  assumed  that  these  collision 
efficiencies  closely  approximate  those  for  thin  hexagonal 
plates  of  ice.  The  hydrodynamic  behaviour  of  an  oblate 
spheroid  was  shown  to  be  essentially  the  same  as  that  of 
a  circular  disk  of  the  same  aspect  ratio  (Pitter  et  al,  1973)* 
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Also,  Jayaweera  and  Cottis  (1969)  and  List  and  Scheraenauer 
(1971)  have  shown  experimentally  that  the  hydrodynamic 
behaviour  of  a  simple  hexagonal  plate  can  be  idealized  with 
sufficient  accuracy  by  that  of  a  circular  disk.  For  these 
reasons,  ice  oblate  spheroids  are  used  in  this  study  to 
approximate  the  simple  hexagonal  ice  plates  in  the 
numerical  simulation  of  the  initiation  of  graupel  growth 
computations.  The  tendency  of  ice  plates  to  accrete 
droplets  at  their  corners  as  proposed  by  Knight  and  Knight 
(1973)  is,  of  course,  not  covered  by  this  type  of 
approximation, 

3,2  The  Mathematical  Model. 


In  order  to  model  the  hydrodynamic  interactions 
between  the  ice  crystals  and  cloud  drops  with  as  much 
physical  realism  and  convenience  as  possible,  the  superpo¬ 
sition  scheme  of  Shafrir  and  Gal-Chan  (1971)  and  Pitter 
and  Pruppacher  (197*+)  is  used.  The  model  neglects  the 
close  boundary  effects,  since  it  assumes  that  each  body 
moves  in  the  stream  caused  by  the  fluid  motion  around  the 
other  body  in  isolation.  The  equations  governing  this 
model  are  described  later,  but  some  physical  aspects  are 
described  as  follows.  This  approach  gives  only  moderately 
accurate  accounts  of  the  hydrodynamic  interaction  of  the 
two  bodies,  especially  when  the  two  bodies  are  very  near 


. 
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to  each  other  or  they  have  a  very  low  relative  velocity. 

However,  the  accelerative  forces  which  occur  during  a 
hydrodynamic  interaction  of  the  falling  bodies  are  small  as 
compared  to  the  product  of  the  particle  mass  and  the  gravitational 
acceleration.  Since  total  energy  associated  with  the  flow 
field  is  approximately  proportional  to  the  mass  of  the  body 
(Steinberger  et  al,  1968),  a  reasonable  parameter  for  evaluating 
the  effect  of  the  flow  field  around  the  smaller  body  on  the 
larger  body  is  the  mass  ratio  of  the  two  bodies.  The  use  of 
the  superposition  method  described  above  would,  therefore, 
be  justified  if  the  mass  ratio  of  the  small  body  over  the 
larger  one  is  small  enough. 

This  argument  also  permits  the  neglecting  o 
the  possibility  of  the  flow  field  around  the  drop  causing 
the  ice  crystal  to  tilt.  The  angular  acceleration  is 
proportional  to  the  torque  and  inversely  proportional  to 
the  moment  of  inertia  of  the  body  which  in  turn  is  propor¬ 
tional  to  the  mass  of  the  body.  The  torque  on  the  spheroid 
is  due  to  the  effect  of  the  flow  field  around  the  sphere 
which  is  proportional  to  the  energy  in  the  wake  of  the 
sphere,  therefore,  the  ability  of  the  water  sphere  to  cause 
the  ice  spheroid  to  tilt  is  approximately  proportional  to  the 
mass  ratio  of  the  bodies.  If  the  mass  ratio  is  always  kept 
small,  significant  tilting  of  the  ice  spheroid  is  unlikely  to 
occur.  This  has  proved  to  be  acceptable  since  the  horizontal 
velocity  of  the  ice  spheroid  calculated  from  the  model  is 
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apparently  zero  for  all  cases  of  study, 

A  thin  oblate  spheroid  falling  freely  in  air 

•would  have  its  broadest  extension  perpendicular  to  the  flow 

so  as  to  maximize  its  drag.  In  such  a  way,  its  axis  of 

symmetry  is  parallel  to  the  axis  of  the  flow.  For  a  body 

dropping  vertically  in  air,  it  is  subject  to  two  forces: 

the  gravitational  and  the  drag  force  due  to  its  motion. 

If  the  air  is  moving,  the  drag  force  on  the  body  would  be 

the  one  due  to  the  relative  motion  of  the  body  and  air. 

— ■  * 

The  force  on  a  body  moving  at  V  and  embedded  in  flow 
moving  at  U*  could  be  approximated  by  -§-CDp*A  |  V-U*|2  in 
the  direction  of  lAv*;  where  A  is  the  projection  area  of 
the  body  perpendicular  to  the  flow,  and  p *  is  the  density 
of  the  medium.  As  the  oblate  spheroid  is  assumed  to  be  in 
steady  state,  the  cross  section  area  of  the  body  normal 
to  the  flow  is  just  a  circle  with  radius  of  the  semi-major 
axis  of  the  spheroid,  i.e.  A  =  Tf  Q  L  •  The  dimensional 
drag  force  Fd  is  then  written  as, 

m  =  i  CpTta^p*  I  v-  uf 

1  I  u  -  v I 

or  Fd  =*£c0nal*p*|  V*-U*|  (V-U*)  . 


The  Reynolds  number  Re-j  is  defined  for  the  above  body  as 


Re 


=  _ L_  ,  so  the  drag  force  can  be 


rewritten  as  Fd  =  -n.  CDReL/)U*ClL*(  v*-  U*) 


Thus,  the 
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trajectories  of  a  large  oblate  spheroidal  ice  crystal  and 
a  small  spherical  water  drop  can  be  described  by  the 
equations  of  motion  as  given  by  Shafrin  and  Gal-Chan  (1971) 
and  Pitter  and  Pruppacher  (197*0: 


dt* 


=  m’5*- JL u*G>LReLa* ( v*  -  U*  ) 

°  4  / 


(3.D 


and 


rr\t 


=  vn*  2*-  2-  Resa  *  (  V*-U  *  ) 


dt* 


(3.2) 


The  first  terms  of  the  right-hand  side  of  the  above 
equations  are  the  gravitational  forces  and  the  second  terms 
are  the  drag  forces  on  the  bodies  due  to  their  relative 
motion  in  air,  N  m  's'  denote  the  masses,'  V*s  1  the  velocities 
j/s  f  the  drag  coefficients  and  Cl^s^the  semi-major  axes 
normal  to  the  flow,  v U*/ s  /  denote  the  flow  fields  due  to 
the  motion  of  the  other  bodies  and  the  subscripts  '  L  '  and 
’ S  '  represent  the  ice  crystal  and  water  drop,  respectively, 
vg* f  and  are  the  gravity  acceleration  and  dynamic 
viscosity  of  air.  The  buoyancy  of  ice  or  water  in  air  has 
been  left  out. 

The  quantities  in  equation  (3.1)  and  (3.2) 
can  be  non-dimensionalized  by  the  relations: 
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v  =  v*/  v‘lM 
u  =  u  V  v‘w 

S=  s*cuVv*i 
t  =  t*v*„/aL* 

/a= 


,  a*  and  V  L00  are  the  dimensional  mass,  the  semi-major 
axis  normal  to  the  flow  and  the  terminal  velocity  of  the  ice 
crystal.  V,  U,  g  and  t  are  the  dimensionless  quantities 
corresponding  to  those  with  asterisks.  During  interaction 
of  the  ice  spheroid  and  water  sphere,  the  Reynolds  numbers 
and  the  drag  coefficients  of  the  particles  vary  due  to 
changes  in  their  terminal  velocity.  The  variation  in  the 
product  ReCj)  for  the  ice  spheroid  is  negligible  since  change 
in  its  terminal  velocity  due  to  the  flow  field  of  the  small 
water  drop  is  negligible.  For  the  small  sphere,  the  drag 
coefficient  times  Reynolds  number  was  found  to  change  by  a 
factor  of  about  "two"  as  it  approached  the  spheroid.  Pitter's 
assumption  that  ReC-p  is  constant  would  be  good  for  Stokes 
flow  where  is  inversely  proportional  to  Re.  However  for 
the  larger  Re  encountered  here  it  appears  to  be  inadequate. 
The  effect  is  more  predominant  if  the  mass  ratio  of  the  small 
water  drop  to  the  large  ice  crystal  is  small.  In  the 
computation  of  the  trajectories,  the  variation  in  drag 
coefficient  and  Reynolds  number  of  the  large  i  e  spheroid 
is  not  be  taken  into  account  but  that  of  the  small  water 
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drop  is  included.  The  correction  is,  therefore,  taken  in 
every  step  of  the  computation  of  the  particle  trajectories 
as  they  approach  the  spheroid.  For  the  large  particle, 
the  drag  coefficient  and  Reynolds  number  are  those  for  a 
body  at  its  terminal  velocity  in  an  undisturbed  medium  and 
for  the  small  particle,  the  drag  coefficient  and  Reynolds 
number  at  its  terminal  velocity  when  falling  in  steady 
state  are  used  to  start  the  computation  of  trajectories  and 
they  are  changed  to  correspond  to  the  variation  in  velocity. 

Applying  the  relations  (3*3),  equation  (3*1) 
can  be  written  in  dimensionless  form  as, 


dV.  =  g  _£  uC0lRcl(Vl'IO  . 


(3.*+) 


Zt  S  4/ 


In  the  case  of  the  ice  spheroid  falling  at  its  terminal 
velocity  in  an  undisturbed  medium: 


d  Vl  =  0  ;  Us  =■  0  and  Vt  -  VL«>  -  1 


dt 


then  equation  (3.*t)  gives 


(3.5) 


For  convenience,  the  parameter  QL  is  equated  to  the  scalar 
magnitude  of  the  non-dimensional  acceleration  of  gravity, 


' 


. 
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i.e. 


Ql  =  i^CotReL  ,  (3.6) 

which  is  dependent  on  variables  which  are  constants  for  a 
given  ice  crystal  and  environment*  Equation  (3*2)  is  also 
non-dimensionalized  with  respect  to  the  terminal  velocity 
and  semi-major  axis  of  large  ice  crystal  for  convenience* 
Since , 


mL=  4-TT<aL!  ArpL/3 

and  1^1^=  4-TTQ.s3ps/3 


where  A/  is  the  aspect  ratio,  then  the  ratio  of  masses  of 
the  two  particles  is 


5  . 

ms  __  as  Ps 
WL  at!ArpL 


Defining  the  ratio  p  =  —  and 

aL 

be  rewritten  as, 

 P 1 


the  mass  ratio  can 


qAr 

Thus,  equation  (3*2)  can  be  expressed  in  non-dimensional 
form  as 


d —  a  ^CpsResq Aril  /  » .  _  | .  \ 
At  5  -  4  pi  L  V$  UL! 


(3.7) 


. 
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For  simplicity,  the  parameter  Qs  is  designated  as 

Os  =  ^  Cos  q  Ar /A  ( 

4  p2- 

where  Qs  is  a  variable  which  depends  on  constants  for  a 
given  ice  crystal,  water  drop  and  atmosphere  and  varies 
with  the  velocity  of  the  water  drop. 


3*3  Formulation  of  Computational  Equations. 


Velocity  is  the  time  derivative  of  position.  In 
the  case  of  either  the  oblate  ice  spheroid  or  the  water  drop, 
the  relation 


(3.9) 


holds  true.  R  is  the  position  vector.  A  system  of  eight 
first  order  differential  equation  can  be  specified  for  the 
water  sphere  and  the  thin  oblate  spheroid  trajectories 
by  resolving  equations  (3*9),  (3.7)-  and  (3#*0  into  their 
components  with  the  simplification  of  equations  (3*6) 
and  (3.8).  Namely,  the  equations  in  their  component 
forms  are: 
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=  V, 


ARli 

dt  VLZ 

#?=  aL[i-(vtz-uSI)] 

dt 


=  V 


LY 


dRLv 

dt 

^=-QL(VLY-Usv) 


dR 


SZ 


dt 


=  v 


SI 


G.-QsCVsZ-Ulx) 

at 

=  “GtS  (  VSv  -  Uly) 

dt 


sr 


*\ 


J 


(3.10) 


Subscripts  2.  and  Y  denote  the  vertical  and  horizontal 
components  respectively*  The  first  four  equations  of  (3*10) 
are  those  for  the  vertical  and  horizontal  time  derivatives 
of  position  and  velocity  of  the  ice  crystal  and  the  last 
four  equations  govern  those  of  the  water  drop. 

The  above  equations  are  solved  numerically 
using  the  Hamming  predictor-corrector-modifier  method* 
(James  et  al,  1967)*  The  Hamming  method  for  numerical 
solution  of  ordinary  differential  equations  is  of  fourth 
order  accuracy*  The  disadvantage  of  this  method  is  that 
it  requires  four  initial  known  step  values  for  calculating 
the  value  of  the  next  step.  To  overcome  this,  the  Runge- 
Kutta  method  (Gerald,  1970)  is  used  to  start  the  integration. 
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The  Runge-Kutta  scheme  is  of  fourth  order  accuracy  also, 
but  it  requires  longer  computational  times.  The  integration 
routine  is  given  in  the  Appendix  III.  During  the  integra¬ 
tion  of  the  trajectories,  it  is  assumed  that  each  body 
experiences  an  external  force  due  to  the  presence  of  the 
other  body  as  though  the  flow  field  of  the  second  body  acts 
at  the  center  of  the  first  body.  The  flow  fields  are  evaluated 
at  the  specific  points  at  each  step  of  the  integration 
by  using  bi-quadratic  interpolation  in  the  appropriate 
co-ordinate  systems.  That  is,  by  taking  sixteen  grid-point 
values  (see  figure  3.1)  around  the  point  in  question.  The 
value  at  that  point  is  found  using  the  expressions  below. 

The  Q’s  are  interpolated  with  the  four  grid-point  values  P’s 
along  the  lines  passing  them  and  the  magnitude  of  X  is  then 
determined  by  the  four  values  of  Q's  with  the  equations, 


Q-  N  -  P*2  ,N  +  "  P2, n)  Sy 


X”  Gl  2  (  Q- s  "  GU )  S  * 

1-  (  GL.+  Gu-  GU-QD(  Sx)/4 


The  flow  fields  due  to  the  motions  of  the  ice  crystal  and 
water  drop  are  calculated  by  using  the  finite  difference 
scheme  as  disscussed  in  Chapter  2  and  the  computer  programs 
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Fig,  3.1 

Illustration  of  Grid-points  for  the  Bi-quadratic 
Interpolation,, 
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are  given  in  Appendix  I  and  Appendix  II  respectively. 

Beyond  the  outer  boundaries  used  to  determine  the  flow  fields, 
the  fluid  is  assumed  to  be  unperturbed  by  the  presence  of 
the  moving  body.  Within  the  boundaries,  the  flow  fields 
are  assumed  to  remain  the  same  as  those  caused  by  the  fluid 
motion  due  to  the  movement  of  the  body  in  isolation  despite 

the  presence  of  the  other  body  -  superposition,  as 

discussed  in  section  3*1  of  this  chapter. 

The  vertical  and  horizontal  seperations  between 
the  centers  of  the  oblate  ice  spheroid  and  the  water  sphere 
are  determined  by: 


Z  =  Rsz  -  Rlz 
Y  =  Rsy  -  R  Ly  , 


(3.H) 


and  these  values  are  assigned  as  the  position  for  evaluating 
the  flow  fields  acting  on  the  other  body  by  using  the 
bi-quadratic  interpolation  method  as  mentioned  above. 


3.^  Evaluation  of  the  Drag  Coefficient,  the  Terminal 

Velocity  and  the  Radius  of  Body  Normal  to  the  Flow 

Le  Clair  et  al  (1970),  after  an  extensive  study 
of  variation  of  drag  coefficient  with  Reynolds  number  of 
spheres,  suggest  the  drag  coefficient  for  a  sphere  of 
Reynolds  number  less  than  0.1  can  be  approximated  by  the 
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expression 


2b 

Re 


(  1 


3 

nr 


Re  ). 


Ever  since  the  study  of  fluid  started,  no  simple  relation 
of  drag  coefficient  with  size  of  an  oblate  spheroid  moving 
in  a  medium  has  been  established.  For  this  reason,  the 
drag  coefficient  of  an  oblate  spheroid  is  calculated  by 
the  method  discussed  as  follows: 

The  dimensionless  local  shear  stress  at  the 
surface  of  a  sphere  X re 

and  Reynolds  number  is  given  by  Masilyah  (1973)  as 


jr=1  in  terms  of  vorticity  ^  r=1 


Tn>U=i 


(3.12) 


2  2 

The  surface  of  the  sphere  for  x=0  is  given  by  Z  +  Y  =  1, 


where 


Z  =  r  cos© 
Y  =  r  sin© 


r=1 


The  surface  area  of  the  circular  strip  in  fig  (3*2)  on 

f 

which  the  shear  stress  is  considered  to  act  is  2  TTvj(dz‘ 
+  dy  f  .  The  drag  force  in  the  direction  of  the 
undisturbed  flow  acting  on  the  whole  surface  due  to  the 
shear  stress  alone  is  then  given  by 


F  =  2  T[\j(dz2  +  dy2  f  T*  e|r-, 


COS  ck 


1 


(3.14) 


. 


K.'*  '  "  n 


■  -•  jj  •  •  r"i  «»W  d9i<B 


* 


51 


Fig.  3.2 

Action  of  the  normal  and  tangential 
stresses  on  an  element  of  area  of  a 
sphere. 


. 
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where  ck  can  be  expressed  as  —  •  The  drag 

on  a  body  is  usually  expressed  in  terms  of  a  dimensionless 
drag  coefficient, 


drag  force  in  the  flow  direction 
(.projection  area  normal  to  the  direction  of 


x  1 _ 

flow)  x  kinetic  pressure* 


Therefore,  the  skin  (friction)  drag  coefficient  can  be 
written  from  equation  (3.1*+-)  as, 


j9  2u  (dzz+dy2)sm0 

TUip’U’2) 


Now  COS  cA.  = 


dy 

(dz^  +  dy^  Y 


5 


dz  =  cos©  d0 
dy  =  sin0  d0 


(3.1 5) 


the  dimension  relation,  x  ="C*/2  p*U*a  ,  and 
making  use  of  equation  (3.12),  equation  (3.15)  yields 


c 


DS 


de 


(3.16) 


The  above  equation  gives  the  skin  drag  coefficient  for  a 
sphere*  For  the  case  of  an  oblate  spheroid,  using  the  same 
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arguments,  equation  (3# 16)  can  be  transformed  to 

fTT 

CDs  =  -A—  tanh^aJ  sin*  ^  ^  d  #  (3.17 

The  pressure  distribution  around  a  sphere  is 
obtained  by  considering  the  radial  and  angular  components 
of  the  Navier-Stokes  equation.  Taking  the  curl  of  the 
equation  and  applying  the  continuity  equation,  the  equation 
(2,4)  becomes: 


(3.18) 


Now  consider  the  angular  component  of  equation  (3.18)  which 
is  given  by: 


At  the  surface,  V*^  =  0,  V*r  =  0  and  for  an  oblate  spheroid 
h.j  =  h£  ,  so  that  equation  (3.19)  becomes 


(3.20) 


Introducing  the  dimensionless  pressure 


P 


P*  -  P*o 


(3.20  a) 


r 
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where  P* 0  is  the  dimensional  frontal  stagnation  pressure, 
substituting  for  h^  and  rendering  all  terms  dimensionless, 
the  equation  (3*20)  becomes: 


(3.21) 


Integration  of  equation  (3.21)  along  the  surface  gives 


+  %  r.-l'j  d  0 


(3.22) 


P  -  Po  = 


where  Po  is  the  dimensionless  frontal  stagnation  pressure 
(at  9=0)  and  P  is  the  dimensionless  pressure  at  any  one 
point  on  the  surface  of  the  sphere. 


The  dimensionless  stagnation  pressure  Po  is 


obtained  by  consideration  of  the  radial  component  of  equation 
(3.18),  Using  the  axisymmetric  properties,  equation  (3.18) 
yields  for  the  radial  direction  as 


Along  the  axis  of  the  spheroid,  V*©  =  0,  hence  equation 
(3.23)  becomes: 


Z±  (<  Cot  0+15]  =  2  w  5>Vr  2B- 


Re  de  e=o  dr  dr 


d  0  0  =  0 


(3.24) 
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■where  all  quantities  have  been  rendered  dimensionless  and 
the  values  of  h^ ,  h^  and  h^  for  a  sphere  have  been 
substituted.  At  the  outer  boundary  P  =  0  and  V=1 
and  at  the  surface  P=Po  and  V=0.  On  the  axis  of  the 
sphere  the  integration  of  equation  (3.24)  from  the  surface 
to  the  outer  boundary,  therefore,  becomes: 


Po  =  1  + 


(3.25) 


Along  the  axis  0=0,  cot  0  =  oo  and  ^  =0,  using  l'Hospital's 
rule,  equation  (3.25)  becomes: 


Po  =  1  + 


Re  Ji  3  9 


dr 

8-0 


(3.26) 


From  figure  (3.1),  the  form  (pressure)  drag  coefficient  is 
given  by: 


C 


DF  = 


$e  2TT  (dz2  +  dy2P  swol  P 


(3.27) 


O  Q  -A- 

where  sine*  =  dz/(dz^  +  dyd  )2  and  at  the  surface,  z=  cos© 
and  dz  =  sin0d0,  rendering  equation  (3*2 7), 


f  TT 

Cdf=  1  P  sih  20  d9 

Jo 


(3.28) 


Equations  (3.22),  (3.26)  and  (3.28)  can  be  transformed  for 
the  case  of  an  oblate  spheroid  as  given  by  Pitter  et  al  (1973) 
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as 


P  -  Po 


Po 


and 


4 


Re 


1  + 


0 


8 


Re 


$ 


11 
?a  s  n. 


-1=0  , 


r* 

?  sin  2 1  dll 


(3.29) 

(3.30) 


(3.3D 


The  total  drag  coefficient  is  obtained  by  adding  the  skin 
and  form  drag  coefficients,  i.e.  Q)  =  Q)S  +  0)F. 


For  a  spherical  body  falling  in  steady  state 
in  air,  the  gravitational  force  acting  on  the  body  should 
be  neutralized  by  the  drag  force  caused  by  the  motion. 
Using  the  definition  of  Re  and  lYls  this  concept  can  be 
written  as 


-~E-CDReasVoosu.  =  w\sg 

t*  / 


CPs -Pa) 

Ps 


(3.32) 


and  can  be  rewritten  as 


(3.33) 


For  known  Reynolds  number  and  drag  coefficient  solving  for 
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the  radius  of  the  sphere,  one  has 

a5  =  (3ya2CDRe2/32g((VpJpa')'*  .  (3.3*0 

The  terminal  velocity  can  then  be  found  using  the  relation 
of  Reynolds  number  to  radius  and  terminal  velocity  of  a 
sphere. 

For  an  oblate  spheroid,  following  the  same 
arguments  with  the  mass,  yy\L  -  -4-  TT  aL3  Ar  pu  ?  thus 

C,Re2  =  32  §  aL5A,(  pL-  pa)  pa  /yu.2  (3<3?) 

or 

aL=  (3  ,U2CDRe752Q  (pL-pa)paAr)^ 

1  5  f  l  (3.35a) 

The  drag  coefficients  for  the  spheroid  and  the 
sphere  in  equations  (3*15)?  (3*28),  (3*17)  are  integrated 
using  finite  difference  method  with  the  vorticities 
calculated  from  the  flow  field  routine.  They  are  then 
used  for  solving  the  radii  of  the  bodies  normal  to  the 
flows  and  the  terminal  velocities  of  the  bodies  which  are 
to  be  applied  to  initiate  the  computation  of  the 
trajectories. 


. 
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CHAPTER  IV  NUMERICAL  PROCEDURES  AND  ASSUMPTIONS, 


4,  1  Numerical  Procedures 

The  flow  field  of  the  ice  crystal  is  evaluated 
with  the  finite  difference  method  discussed  in  section 
(2,5.1)  of  Chapter  2,  The  grid  points  used  are  in  oblate 
spheroidal  co-ordinates  with  radial  step  size  A^=  A  =0.1 


and  angular  step  size 


Boundary  conditions 


shown  in  section  (2,5*3)  are  used.  Outer  boundaries 
are  located  at  50  crystal  semi-major  axis  lengths. 
Approximate  grid-point  values  are  evaluated  from  the  analy¬ 
tical  solutions  of  the  Stokes  flow  around  a  sphere 
(equation  2,14),  A  relaxation  technique  is  then  employed 
to  acheive  more  accurate  grid  values.  The  solutions  are 
considered  to  have  converged  when  successive  values  at  all 
grid  points  for  both  the  dimensionless  stream  function 
and  the  non-dimensional  vorticity  fields  have  fractional 
changes  of  less  than  0,001,  An  over-relaxation  method 
is  used  for  equation  (2,39)  while  an  under-relaxation  is 
applied  to  equation  (2,38),  as  the  former  equation  is 
mathematically  linear  and  the  second  is  not. 


Two  methods  are  used  to  calculate  the  flow 


field  of  the  water  drop.  For  droplets  of  Reynolds  number 
greater  than  or  equal  to  0,1,  the  finite  difference  method 
of  solving  the  complete  Navier-Stokes  equation  numerically 
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as  discussed  in  section  (2*5*2)  is  used  at  grid-points 
in  modified  spherical  co-ordinates  with  radial  step  size 
C  =  0.05  and  angular  step  size  D  =  3°  •  Boundary  conditions 
shown  in  section  (2.5**+)  are  used  with  outer  boundaries 
located  at  80  drop  radii*  Approximate  initial  grid-values 
are  taken  from  analytic  solution  of  creeping  flow  past  a 
spherical  obstacle  (equation  2*14),  Relaxation  techniques 
and  convergence  criteria  are  followed  the  same  as  that  for 
the  spheroid.  For  drops  of  Reynolds  number  smaller  than 
0.1,  the  flow  fields  are  computed  using  a  procedure  suggested 
by  Le  Clair  et  al  (1970).  This  involves  the  Stokes  solution 
for  a  bounded  fluid. 

The  flow  fields  evaluated  from  the  above 
discussions  are  in  terms  of  dimensionless  stream  functions 
and  vorticities.  These  values  have  to  be  converted  into 
non-dimensional  velocities.  In  two-dimensional  curvilinear 
co-ordinates,  the  dimensional  velocities  can  be  related  to 
the  dimensional  stream  function  by  the  expressions: 

v*  =  -Ml 

K  2  ^5 

and  v*  =  _L_  1^! 
q2  h,hj  dq, 

h^  ,  h^  ,  h^  are  defined  as  the  same  in  Chapter  2.  Therefore, 
in  spheroidal  co-ordinates,  the  non-dimensional  velocities 
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can  be  expressed  as: 


V, 


-  cos 


and 


V, 


_ _ _ 

(s»'nh2^  + cos2^)cosh^$in’^ 

_ cosk2  _  3  ^ 

(sinh2£;  +•  cos2  ;|)  cosK^  sen  V|  3  ^ J  ^ 


>(4.2) 


and  in  the  modified  spherical  co-ordinates  with  r=e  , 
dimensionless  velocities  are 


the 


and 


-e‘2Z 

3^ 

sin  Q 

ae 

e-£Z 

sin  0 

az  J 

(*+•3) 


In  order  to  save  computing  time,  the  symmetry  condition 
around  the  OL-axis  is  used,  thus  implying  that 


and 


r[=o(Ti 


~  ^ 


V 


-o.n 


3Ve 
3  0" 


6=-0,7T  -  0  . 


(4.4) 


(4.5) 


Using  1 'Hospital's rule  and  equations  (4.2)  and  (4.3), 
ye  have: 


Vq  1 0, tj 

V  ^  I  0, 11 


=  0 

g  -  cosh2Ca  32^ 
cosh1  ^  3  Kl1 


(4.6) 
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and 


le-o,Ti 

Vr  |  e  =o,  r 


d  e1 


* 


(4.7) 


The  flow  fields  calculated  from  the  previous 
arguments  are  for  stationary  obstacles  in  a  moving  stream. 
Por  a  body  moving  in  a  calm  fluid  the  velocities  would  be 
relative  to  those  calculates  above.  Thus,  the  velocities 
would  be  the  differences  between  the  terminal  velocity  of 
the  body  and  those  calculated  from  the  case  of  a  stationary 
obstacle  in  a  moving  medium. 

The  corrected  flow  fields  are  then  applied 
into  the  trajectory  calculations.  In  the  trajectories 
computation  routine,  the  ice  spheroid  and  the  water  drop 
are  initially  vertically  separated  by  three  semi-major  axes 
of  the  crystal.  This  separation  is  believed  to  be  enough 
since  the  change  in  the  velocity  of  the  flow  due  to  the 
spheroid *s  approach  is  practically  negligible  at  a  distance 
of  three  semi-major  axis  upstream.  The  horizontal  separation 
of  the  fall  paths  of  the  two  particles  is  varied  for  each 
calculation.  Therefore,  the  initial  horizontal  and  vertical 
acceleration  of  the  ice  crystal  and  the  water  drop  are 
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taken  to  be  zero  and  the  vertical  velocities  are  their 
terminal  velocities  non-dimensionalized  by  the  terminal 
velocity  of  the  ice  crystal.  Initially,  the  flow  field 
caused  by  one  body  at  the  location  of  the  other  body  is 
taken  to  be  zero.  The  bodies  are  initially  falling  vertically 
at  their  terminal  velocities.  That  is,  the  horizontal  velocities 
are  taken  to  be  zero  at  the  beginning  of  the  integration. 

The  non-dimensional  time  difference  used 
between  every  step  in  the  integration  of  the  trajectories 
is  0.01.  For  an  ice  crystal  of  semi -major  axis  of  400 
microns  and  falling  at  a  speed  of  45  cm  sec~  ,  the  actual 
time  step  would  be  less  than  0.1  msec.  This  amount  of 
time  has  proved  to  be  small  enough  to  avoid  errors  due  to 
the  numerical  integration. 

The  hydrodynamic  interaction  process  is  assumed 

o 

to  be  taking  place  in  an  atmosphere  of  temperature  -IOC  and 

of  pressure  700  mb.  The  density  of  air  is  therefore 

-4 

9.267  x  10  gm  cm.  and  the  dynamic  viscosity  of  air  is 

-4 

I.667  x  10  poise.  The  density  of  ice  is  taken  as  0.92 
gm  cm  and  the  density  of  water  is  1  gm  cm  J  .  After 
proper  non-dimensionalization,  these  values  are  used  with 
appropriate  Reynolds  number,  radii  of  particles,  drag 
coefficients  and  terminal  velocities  to  obtain  the 
parameters  for  the  computation  routine. 
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*+•2  Assumptions  of  the  Numerical  Simulation  of  the 
Initiation  of  Graupel  Growth. 

There  is  no  perfect  modeling  that  can  fully 
illustrate  the  happening s  in  a  graupel  growing  cloud.  In 
every  cloud  model  or  graupel  or  hail  growth  model,  there 
must  be  some  limitations  and  assumptions  involved.  There 
is  no  exception  in  this  model.  As  it  is  impossible  to 
include  a  random  sample  of  size  distribution  of  cloud 
droplets  and  ice  crystals  in  the  model,  the  model  is 
restricted  to  a  particular  ice  crystal  and  it  assumes  that 
drops  of  the  same  radius  are  distributed  everywhere  in 
the  region  under  consideration.  This  assumption  implies 
that  no  collision  between  drops  occurs.  If  collision 
between  drops  exists,  then  the  size  of  drops  can  no  longer 
be  the  designed  one.  Moreover,  that  water  drop  must 
retain  its  spherical  shape  during  the  whole  process  of 
calculation,  along  its  trajectory  and  after  collision  with 
the  ice  crystal.  In  another  words,  the  water  drop  must 
be  spherical  all  the  time  and  it  must  freeze  completely 
without  any  distortion  in  shape  when  it  touches  the  surface 
of  a  crystal.  The  drop  will  remain  at  the  point  of  impact. 
This  is  not  likely  to  happen  exactly  in  the  atmosphere, 
because  according  to  Maruyama  (1968)  the  frozen  cloud 
droplets  are  found  to  have  accreted  together  with  one  another 
so  tightly  that  no  vacant  gaps  are  seen  at  the  points  of 
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adherence  in  the  interior  structure.  This  suggests  that 
■water  drops  do  not  freeze  completely  right  after  the  impact, 
but  some  of  the  water  remains  at  its  surface  until  another 
cloud  droplets  accrete  on  them.  Besides,  results  of 
density  measurements  suggest  that  cloud  droplet  accretion 
is  nearly  uniform  during  the  growth  from  a  small  graupel 
to  a  large  one.  According  to  Mason  (I960),  the  density  of 
a  graupel  is  around  0,6  gm  cm”3  0r  less.  This  density 
is  much  smaller  than  that  of  pure  ice.  This  implies  that 
air  bubbles  should  be  trapped  inside  a  graupel. 

As  particular  interest  is  on  critical  collision 
parameters,  only  these  are  calculated.  That  is,  the 
range  of  discrepancy  between  the  lines  of  motion  of  the 
crystal  (oblate  spheroid)  and  the  sphere  which  results  in  a 
collision  is  evaluated.  The  outer  (or  inner)  most  critical 
point  of  impact  is  defined  as  that  critical  offset  on  the 
ice  crystal  that  away  from  (or  inside  of)  which  collision 
never  occurred  and  inside  of  (or  away  from)  which  collision 
always  occurred  for  particular  size  of  cloud  droplets. 

After  a  critical  point  of  impact,  either  inner  or  outer, 
has  been  calculated,  the  critical  point  on  which  a 
second  drop  would  just  touch  the  first  drop  is  evaluated. 

It  is  assumed  that  the  flow  field  around  the  ice  crystal 
does  not  change  even  when  there  is  a  water  drop  frozen  on 
its  surface.  This  is  the  most  agreeable  assumption.  If 
this  change  in  flow  field  were  taken  into  account,  the  new 
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flow  field  would  have  to  be  calculated  again  from  the 
complete  Navier-Stokes  equation.  Then  the  inner  boundaries 
would  no  longer  be  the  shape  of  a  simple  oblate  spheroid. 

It  becomes  a  combination  of  an  oblate  spheroid  and  a  small 
sphere.  The  oblate  spheroidal  co-ordinate  system  would  no 
longer  provide  the  convenience  it  gave  before.  No  simple 
co-ordinate  system  is  suitable  for  an  obstacle  of  this 
shape.  Moreover,  the  computer  time  used  to  generate  a 
flow  field  around  a  symmetric  obstacle,  by  solving  the 
Navier-Stokes  equation  numerically  using  the  finite  difference 
method,  is  very  long  already.  It  would  take  even  much  longer 
to  obtain  the  flow  field  around  an  irregular  particle.  To 
grow  a  whole  graupel  is  not  the  object  of  this  study. 

The  object  of  this  model  is  to  determine  the  maximum  angle 
that  accreted  water  drops  would  form  with  the  ice  crystal. 

Two  or  three  drops  adhered  at  the  outer  most  or  inner  most 
side  of  the  ice  crystal  would  be  good  enough  to  illustrate 
the  solution.  If  a  few  small  drops  are  attached  to  the 
surface  of  an  ice  crystal,  the  flow  field  around  it  would 
not  be  changed  a  very  great  deal  from  that  of  a  clean  and 
smooth  one  of  the  same  size.  For  all  these  reasons,  the 
flow  field  around  the  water  drop  adhered  ice  crystal  is 
taken  from  that  of  an  ice  oblate  spheroid  of  the  same 
Reynolds  number  and  aspect  ratio. 

It  has  been  mentioned  in  the  Chapter  3  that 
the  drag  coefficient  and  Reynolds  number  of  the  ice  crystal 
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are  assumed  to  be  constant  throughout  the  integration  of 

the  trajectories  of  the  two  particles  though  there  is  a 

n 

very  small  deceleration  of  order  of  10  are  reported 
from  the  model.  In  addition  to  all  the  assumptions  discus¬ 
sed  above,  the  buoyancy  effect  of  ice  and  water  in  air  are 
neglected  in  the  calculations. 


.  ,  -  . 


CHAPTER  V  NUMERICAL  RESULTS  AND  DISCUSSION 


5. 1  The  Flow  Fields  around  an  Oblate  Spheroid  (Ice  Crystal). 

In  most  continental  clouds,  the  drops  are  too 
small  to  grow  by  collision  and  coalescence  with  other 
drops.  In  such  clouds,  when  portions  of  cloud  become  cooled 
enough,  ice  crystals  are  nucleated.  The  ice  crystals  formed 
usually  grow  in  two  stages:  initially  by  diffusion  of  water 
vapour,  and  subsequently  by  riming.  Ono  (1969)  and  Wilkins 
and  Auer  (1970)  showed  that,  after  nucleation,  thin 
hexagonal  ice  plates  must  grow  by  vapour  diffusion  to  at 
least  200  microns  in  radius  before  riming  may  commence. 

For  this  reason  and  due  to  the  limitation  of  time  and  comput¬ 
ing  funds,  two  sizes  of  ice  oblate  spheroids  are  studied, 
namely  with  Reynolds  number  10  and  20,  which  correspond 
to  radii  of  289.2  and  39603  microns.  The  drag  coefficients 
obtained  from  the  converged  numerical  solution  of  the 
Navier-Stokes  equation  of  motion  for  viscous  flow  past  an 
oblate  spheroid  of  axis  ratio  0o05  of  Reynolds  number  10 
and  20  are  2.81  and  3*92.  The  stream  function  and  vorticity 
fields  near  the  spheroid  are  plotted  for  Reynolds  number  20 
in  figures  5.1.1  and  5.1.2,  and  for  Reynolds  number  10  in 
figures  5.1.3  and  5.1.^.  The  computed  values  of  the 
vorticity  at  the  surface  of  the  body  of  Reynolds  number  20 
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Pig.  J.1.1 

Dimensionless  Stream  function  of  an 
of  Re  =  20  and  Ar  =  0.05 


3-0 


Oblate  Spheroid 
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Pig.  5«1« 2 

Dimensionless  Vorticity  Field  Around  An  Oblate 
Spheroid  of  Ar  =  0.05,  Re  =  20 


70 


Fig.  5#l« 3 

Dimensionless  Stream  Function  Field  Around  An 
Oblate  Spheroid  of  Ar  =  0.05  Re  =  10 
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Pig.  5.1. ^ 

Dimensionless  Vorticity  of  Oblate  Spheroid  of 
Re  =  10  ,  Ar  =0-05 
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and  10  together  with  those  from  Pitter  et  al  (1973)  for 
Ar  =  0.05  and  those  from  Masilyah  and  Epstein  for  Ar  =  0.2 
are  shown  in  figure  5* 1*5*  After  converting  from  the  stream 
functions,  the  vertical  and  horizontal  velocity  fields  of 
air  passing  the  oblate  spheroid  of  Reynolds  number  10 
relative  to  free  air  due  to  the  motion  of  the  particle  are 
plotted  in  figures  5,1,6  and  5*1*7*  The  vorticity  on  the 
surface  of  the  oblate  spheroids  seem  to  agree  pretty  well 
with  those  of  Pitter  et  al  and  the  dimensionless  stream 
functions  at  the  windward  side  of  the  oblate  spheroid  of 
Re  =  20  is  almost  the  same  as  those  from  their  result. 
However,  there  are  some  small  deviations  in  the  stream 
function  field  at  the  downstream  side  of  the  oblate  spheroid. 
This  deviation  is  most  probably  due  to  the  difference  in 
the  converging  criteria  and  the  computation  routine. 

Anyway,  a  very  well  developed  eddy  at  the  downstream  end 
of  the  oblate  spheroid  both  for  Re  =  20  and  Re  =  10 
can  be  seen.  The  eddy  length  —  the  ratio  of  the  length 
of  the  wake  to  the  major  axis  of  the  oblate  spheroid  — 
of  the  present  results  are  plotted  in  figure  5* 1«8  in 
conjunction  with  the  curve  given  by  Pitter  et  al  for  oblate 
spheroids  of  Ar  =  0.05  and  various  Reynolds  numbers. 

It  is  clear  from  the  diagram  that  the  eddy  length  of 
the  present  results  agree  quite  well  with  those  obtained 
by  the  previous  workers. 
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Pig,  5*^*5 

Variation  of  the  Surface  Vorticity  -with 
Polar  Angle  and  Reynolds  Numbers  for 
Oblate  Spheroids 
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Pig.  5*1,6 

Dimensionless  Vertical  Velocity  of  Air  Passing  Oblate 
Spheroid  of  Re=10. 
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Fig.  5.1.7 

Dimensionless  Horizontal  Velocity  of  Air 
Passing  an  Oblate  Spheroid  of  Re  =  10 
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Fig,  5* 1*8 

Variation  with  Reynolds  Number  of  the  Length  of 
the  Standing  Eddy  at  the  Downstream  End  of 
Oblate  Spheroids  of  Axis  Ratio  0,05 
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5. 2  Flov  Fields  around  Spherical  Drops. 


Wilkins  and  Auer  (1970)  in  their  field  study  of 
the  riming  properties  of  hexagonal  ice  crystals  showed  that 
the  size  distribution  of  droplets  accreted  on  ice  crystals 
of  radius  smaller  than  *+00  microns  lies  in  the  range  of 
5  to  18  microns  in  radius.  Fletcher  (1969)  suggests  that 
the  droplets  in  continental  convective  clouds  have  radii 
mainly  smaller  than  10  microns  in  range.  Therefore,  only 
spherical  drops  of  Reynolds  numbers  smaller  than  0.1 
(radius  less  than  17  microns)  are  considered  in  this  study. 
Figures  5*2.1  and  5*2.2  represent  the  dimensionless  stream 
function  and  vorticity  fields  of  air  passing  a  sphere  of 
Reynolds  number  of  0.1.  A  summary  of  the  radii,  drag 
coefficients  and  terminal  velocities  of  spheres  of  various 
Reynolds  numbers  in  an  atmosphere  of  -10°C  and  700  mb  are 
given  in  Table  5*2.1. 


Reynolds 

Number 

10^  x  radius 
(cm) 

Terminal 
Vel  (cm 
sec” ' ) 

Drag  coeff 
from  numer- 
cal  method 

Drag  coeff. 
Le  Clair ’ s 
approx. 

0.1 

19.14 

4.70 

243.50 

244.5 

0.01 

8.83 

1.02 

2399*93 

2404. 5 

0.005 

7.01 

0.64 

4802.89 

4804.5 
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0.004  6.51  0.55  6010.93  6004.5 

0.003  5.91  0.40  8001.38  8004.5 

Table  5.2.1  Reynolds  Number,  Drag 
Coefficient,  Radius  and  Terminal 
Velocity  of  Spherical  VJater  Drops 
in  Atmosphere  of  j-10  C,  7Q0  mb 
(  =  9.267  x  lor4'  gm  cm-^, 

=  1 -667  X  10“  poise  ) 

Table  5*2.1  shows  the  two  drag  coefficients  calculated  by 
the  method  discussed  in  section  3*3  using  the  analytic 
flow  field  solution  for  Stokes  flow  and  the  relation  of 
drag  coefficient  and  Reynolds  number  suggested  by  le  Clair 
et  al  (1970)  for  small  speres  of  Reynolds  number  less  than 
0.1.  The  differences  between  the  two  methods  are  well 
within  1  percent. 

5. 3  Trajectories  of  Water  Drops  Relative  to  an  Ice 
Oblate  Spheroid. 

The  superposition  method  is  employed  with  the 
ice  spheroids  of  Reynolds  numbers  of  10  and  20,  and  water 
spheres  of  Reynolds  number  of  0.003,  0.004,  0.005j  0.01  and 
0.1.  The  trajectories  of  a  water  sphere  of  5*91  microns 
and  of  1 9 • 1 4  microns  radius  (Reynolds  number  0.003  and 
0.1  respectively)  interacting  with  an  ice  spheroid  with 
semi-major  axis  of  396.3  microns  (Reynolds  number  20)  are 
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Pig,  5.2.1 

Dimensionless  Stream  Function  of  Sphere  of  Re  =  0.1 


.  . *•  •  • 


-■ 


o.oi 


Pig,  5.2.2 

Dimensionless  Vorticity  of  a  Sphere  of  Re  =  0.1 


'  •  -  *  .  . 


’ 


. 


81 


plotted  in  figures  5«3*1  and  5*3* 2.  The  curvature  of  the 
trajectories  in  all  the  cases  of  this  study  is  not  as  great 
as  those  obtained  by  Pitter  et  al  (197*+)  and  the  annular 
behaviour  —  no  accretion  of  droplets  in  the  center  of 
the  ice  crystal,  even  for  drops  as  small  as  5*91  microns 
in  radius,  is  not  so  predominant*.  This  would  be  mainly  due 
to  the  correction  of  the  parameter  Qs,  in  equation  3*10, 
as  a  result  of  the  change  of  velocity  of  the  water  drop* 

When  the  ice  oblate  spheroid  approaches  the  water  drop,  the 
velocity  of  the  drop  increases  due  to  the  fluid  motion 
immediately  surrounding  the  ice  particle.  As  a  result,  the 
Reynolds  number  of  the  water  sphere  increases  proportionally. 
The  Best  number  CoRe  of  a  body  moving  in  a  fluid  is  a 
constant.  If  the  Reynolds  number  of  the  Water  drop  increases, 
the  product  of  the  drag  coefficient  and  the  Reynolds  number 
decreases,  i«e.  the  parameter  Qs  reduces  in  magnitude. 

Then  the  terminal  velocity  of  the  water  sphere  is  never 
able  to  increase  to  the  terminal  velocity  of  the  oblate 
spheroid  before  collision.  Therefore  the  horizontal  forces 
are  not  great  enough  to  move  a  water  sphere  which  approaches 
the  center  of  the  spheroid  around  the  oblate  spheroid. 

In  Pitter  et  al*s  study,  the  Reynolds  number  of  the  water 
sphere  is  assumed  constant  throughout  the  whole  interaction. 
The  parameter  Qs  is,  therefore,  held  constant.  The  vertical 
and  horizontal  acceration  of  the  water  sphere  is  then 
capable  of  becoming  large  enough  to  enable  the  drop  to  speed 
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Pig.  5.3.1 

Trajectories  of  5.91  micron  Water  Drop  Relative 


to  an  Oblate  Spheroid  of  ice  of  Re  =  20 
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Trajectories  of  19*1*+  microns  Water  Drop 
Relative  to  an  Oblate  Spheroid  of  Ice  of  Re  =  20 
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up  and  to  go  around  the  oblate  spheroid  of  ice  without 
collision. 

The  maximum  angle  that  water  drops  of  radius 
5.91  microns  formed  when  they  collected  at  the  outer  edge 
of  an  oblate  ice  spheroid  of  Reynolds  number  20  is  drawn 
in  figure  5*3*3*  Figure  5’*3*1+  represents  the  variation 
of  the  outer  maximum  angles  that  water  drops  may  create 
on  oblate  spheroids  of  Reynolds  number  10  and  20  with 
various  Reynolds  numbers  of  water  spheres.  Since  no  water 
droplet  which  approaches  the  middle  part  of  the  ice  crystal 
expelled  around  without  collision  occured,  the  inner 
critical  angle  would  no  longer  be  of  any  importance  and 
they  are  no't  recorded.  The  outer  critical  angle  that  water 
sphere  would  create  on  an  oblate  ice  spheroid  is  rather 
large,  but  they  are  the  maximum  angle  that  the  drop  just 
touched  the  other  drop.  However,  in  the  atmosphere,  the 
second  drop  most  likely  cannot  attach  to  a  single  drop 
adhered  to  the  oblate  spheroid  of  ice  by  just  touching  it. 
Therefore,  the  vertex  angle  that  a  graupel  would  have 
should  be  smaller  than  twice  the  maximum  angle  given  above. 
The  bigger  the  ice  crystal  and  the  smaller  the  cloud  water 
droplet,  the  sharper  would  be  the  angle  of  the  point  of  the 
graupel. 

In  all  the  cases  of  this  study,  no  capture  of 
water  drops  in  the  wake  of  an  ice  oblate  spheroid  is  ever 
detected.  Examining  the  vertical  and  horizontal  isotachs 
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Drops  Forms  on  an  Oblate  Spheroid  of  Re 
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plots  in  figures  5*1*6  and  5*1*7,  when  a  very  small  water 
sphere  passes  the  ice  spheroid,  there  is  a  tendency  to  move 
the  water  sphere  towards  the  axis  of  the  ice  particle  and 
the  vertical  flov;  field  at  the  rear  of  the  ice  body  may 
push  the  water  drop  downward  (i.e  .  back  toward  the 
trailing  end  of  the  spheroid)*  If  the  parameter  Qs  in 
equation  3*10  still  remains  comparably  very  large  after 
passing  around  the  ice  spheroid,  it  may  be  able  to  accelerate 
to  a  velocity  to  catch  up  with  the  larger  particle  which 
is  in  fact  momentarily  decelerating  very  slowly  due  to  the 
flow  field  of  the  water  drop*  It  will,  however,  soon  regain  its 
terminal  velocity  again.  This  may  explain  why  the  tip  of  a 
graupel  is  seen  in  laboratory  studies  to  be  formed  by  very 
small  cloud  droplets. 
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CHAPTER  VI  CONCLUSIONS  AND  RECOMMENDATION 

The  model  made  in  this  study  was  fairly 
successful  though  it  only  dealt  with  two  dimensionals. 
Convergence  to  the  steady-state  Navier-Stokes  equation 
was  slow.  The  computing  time  was  therefore  very  lomg, 
up  to  30  minutes  CPU  time  in  the  IBM  360  computer  of 
the  University  of  Alberta,  1581  grid-point  values  for 
the  thin  oblate  spheroids  of  ice  and  4941  grid-point 
values  for  the  spherical  water  drop  were  eveluated  to 
get  the  flow  field.  However,  the  model  preliminarily 
revealed  the  following: 

(1)  The  length  of  the  wake  formed  at  the  down-stream 
end  of  the  ice  crystal  is  longer  for  higher  Reynolds  numbers* 

(2)  For  small  Reynolds  number  (5.1  ),  the  drag 

coefficient  can  be  successfully  approximated  by  the 
expression:  Cd~  (  i  . 

(3)  Water  drops  of  Reynolds  number  0*003  to  0*1 
approaching  a  thin  oblate  spheroid  of  ice  of  Reynolds  number 
10  or  20  would  not  pass  around  the  ice  particle  without 
collision* 

(4)  This  study  suggests  that  ice  crystals  appear  to 
accrete  droplets  primarily  on  their  upstream  side.  This 
supports  the  theory  of  graupel  growth  suggested  by  List 
as  discussed  in  section  1*3* 
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(5)  No  wake  capture  was  found  in  this  study, 

(6)  The  maximum  angle  that  water  drops  would  accrete 
on  an  ice  crystal's  edge  depends  on  both  the  size  of  the 
crystal  and  that  of  the  droplet.  Small  crystals  and  big 
droplets  give  larger  angles  of  riming  growth  beyond  the 
edge  of  the  ice  crystals. 

Since  only  mininal  annular  behaviour  was 
detected  when  an  ice  crystal  was  approximated  by  a  thin 
oblate  spheroid,  the  theory  of  the  formation  of  conical 
graupel  at  the  points  and  the  edge  of  hexagonal  plates 
as  hypothesized  by  Knight  and  Knight  could  not  be  further 
justified.  The  collision  efficiency  is  no  doubt  greater 
along  the  edge.  Trajectories  showing  droplets  passing 
around  the  ice  particle  as  obtained  by  Pitter  were  not 

found.  The  reason  of  this  discrepancy  was  explained  earlier. 

In  order  to  fully  understand  the  growth  of 

graupel,  a  three  dimensional  model  should  be  used  instead  — 
that  is,  a  hexagonal  plate  should  be  used  rather  than  an 
oblate  spheroid.  The  boundary  conditions  at  the  corners  and 
along  the  sharp  edges  of  the  hexagonal  plate  may  first 
have  to  be  investigated  by  experimental  technique.  This 
is  important  as  the  corner  effects  would  affect  the  flow 
field  drastically.  Another  numerical  method  should  be 
used  instead  of  the  finite  differencing  method  used  in  this 
study.  The  finite  differencing  scheme  is  rather  slow  and 
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a  convenient  co-ordinate  system  is  hard  to  seek  for  a 
hexagonal  plate.  The  best  alternate  numerical  method  would 
be  the  finite  element  method.  This  scheme  deals  in 
infinitesimal  area  but  not  in  grid-points.  It  is  faster  in 
converging  and  it  can  deal  with  a  three  dimensional  object 
of  any  shape  in  simple  Cartesian  coordinate  system.  The 
superposition  model  would  be  still  valid  provided  the  mass 
ratio  is  small.  The  flow  field  around  the  ice  crystal, 
after  a  water  drop  is  adhered  on  it,  would  be  better 
calculated  again.  For  the  finite  element  method,  this 
procedure  would  not  be  a  hazard  since  the  scheme,  as 
mentioned  before,  is  suitable  for  an  object  of  any  shape 
and  the  computation  time  is  far  shorter  than  the  finite 
differencing  method.  If  possible,  the  range  of  water 
droplets  under  study  should  be  enlarged  so  that  the  feasibi¬ 
lity  of  wake  capture  would  be  fully  investigated. 


■ 
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C  CALCULATION  OF  THE  FLOW  FIELD  AFCUND  AN  OBLATE  SPHEROID 

DIMENSION  Cl  (61)  , C2  (6 1 ) , S INHA  (61)  , COSHA (61)  , 

* SINH2  (51)  ,C3  (31)  ,C4  (31 ) ,SINA  (40)  , COS A  (4  0)  , 

*  COS  2  A  (3  1  )  ,  FEL  (31,51)  ,  GL  (3  1  r  51 )  ,  W  (31,51)  , 

*VO  (31,51)  , SIN 2A  (31)  ,COSH2  (61)  ,KACRC  (61)  , RACES  (31) 

*  ,  ?2  (40)  ,  PT  (40) 

RE=10. 

REFF=1 . 3 

R  E  F  I  =  0 . 001 
RSB=0. 001 
A  =  0. 1 

B  =  6.*3.  1416/180. 

AR=0 . 05 

ZI ?  =  A  LOG  (S  Q?  T  (  ( 1  .  +  A  R)  /  (1 .  ~A  R)  )  ) 

S  IN  HZ  =  S  IN  H  ( Z I A ) 

COSHZ=COSH  ( Z I A ) 

S  ECHZ  =  1 . /COS  HZ 
SECH2  =  SECHZ*SECI1Z 
COSHS=COSHZ*COSHZ 
SINHS  =  S  IN  HZ5*  SIN  HZ 
C5=2. /A/A+2. /B/B 
RAC=?-/C?/3. /COSUZ/A/B 
DO  1  J  = 1  ,51 
ZI=ZIA+  FLOAT  (J- 1 )  * A 
SINHA  (J) =  S I N H  (ZI) 

COS HA  (J) =CCSH  (ZI) 

COSH2  (J)  =SIN HA  (J)  *SINHA  (J) 

SINH2  (J)  =5 INHA  (J)  *SINHA  (J)  • 

Cl  (J)  =  ( 2  .  -  A  *  S I N  H  A  (J)  /COSHA  (J)  )  /2  .  /A/C5/A 
C2  (J)  =  (2.  +A*SINHA  (J)  /COSHA  (J)  )  /2  .  /A/A/C5 
RACRC  (J)  =  ? A C /C 0 S H A  (J) 

1  CONTINUE 

DO  2  1=1,31 
3NU=3*FL0AT  (I- 1 ) 

SINA (I)  =  S I N (3NU) 

COSA  (I)  =COS  (3NU) 

S  IN  2  A  (I)  =S  IN  A  (I)  *5  IN  A  (I) 

COS2A  (I)  =  COSA  (I)  *  C  0  S  A  ( I) 

2  CONTINUE 

DO  10  1=2,30 

C3  (I)  =  (2.  -  3/SIN A  (I)  *C0  5A  (I)  )  /2.  /B/3/C5 
C 4  (I)  =  (2.  +3/SINA  (I)  *COSA  (I)  )  /2.  /B/B/C5 
RACES  (I)  = RAC /SIN A  (I) 

10  CONTINUE 

DO  3  1=1,31 

READ  (1 , 999)  (EEL  (I, J)  , J  =  1 , 51) 

3  CONTINUE 

DO  4  1= 1 , 31 

READ  (2 , 999)  ( VO  (T  ,  J )  ,  J=  1 , 5  1 ) 

4  CONTINUE 

DO  5  J= 1 , 51 
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DO  5  1=1,31 
GL  (I  ,  J)  =V0  (I  ,  J)  *COSHA  ( J)  *  S I N  A  (I)  /CCSH2 
5  CONTINUE 

I FLAG=1 
DO  6  N= 1,500 
L=0 

IF  (1- IF LAG)  91,90,91 
90  DO  8  1=2,30 

IF  (1-2)  30,31,30 

31  DO  21  J  =  2 , 50 

?FEL=F£L(I,J+1)  *C  1  (J)  +  FEL  (I,  J-1)  *C2  (J)  +  FEL  (I  +  1,J)*C3(I) 

*  +  FEL  (1-1 , J) *CU  (I) -GL  (I,J) *SECH2/C5*  (SINH2 (J) +COS2A (I)  ) 

D 1L=RFEL-FZL  (I, J) 

FEL  (I, J) =  FEL  (I, J) +D1 L*33FF 

3 GL  =  GL  (I,  J  +  1 )  *C1  (J)  +GL  (I,  J-1 )  *C2  (0)  +  G  I.  ( I  +  1  ,  J )  *  C  3  (I) 

*  +  GL  (1-1  ,  J)  *C4  (I)  -  RACRC  (J)  *SINA  (I)  *  (  (FEL  (I,  J+1)  -FEL  (I,  J-1)  ) 

**  (GL  (1+  1  ,J)  /SIM2A  (1+1)  -GL  (I,  J)  /3/SIN A  (I)  )  )  +EACRS  (I)  *COSHA  (J) 
**  (  (FEL  (1  +  1 ,  J)  -  FEL  (1-1  ,  J)  )  *  (GL  (I,  J  +  1 )  /COSK2  (J  +  1)  -GL  (I,  J-1) 
*/C0S  H2 (j-1)  )  ) 

D2 L  =  KGL - GL  (I , J) 

GL  (I,  J)  =GL  (I ,  J)  +D2L*RSFG 

IF  (ABS  (  D  1  L)  .  LT.  .  001 .  AND.  ABS  ( D2L)  .  LT .  .  00  1 )  L  =  L  + 1 


21 

CON 

TINUE 

GO 

TO  3 

3  0 

-  r 

(1-30) 

32,33,32 

32 

DO 

22  J  =  2 

,50 

RFL 

L=FE  L  ( 

I, J+1) *C1  (J)  +  FEL  (I,  J-1)  *C2 

(J) 

+  FEL  (I 

+  1 , J)*J3  (I) 

*  +  FE 

L  (1-1 , 

J) ,*C4  (I)  -GL  (I ,  J)  *SECH2/C5* 

(31 

Nil 2  (J) 

+  C052A  (I)  ) 

22 


2  3 
P 


J  *L=EFEL-FEL  (I, J) 

FEL  (I, J) =FE1  (I,  J)  +D1 L*  FEF  F 

?GL=GL  (I,  J  +  1 )  *C1  (J)  +  GL  (I ,  J-  1  )  *C2  (J)  +GL  (1+  1  ,  J)  *0  3  (I) 

*  +  GL  (1-1  ,  J)  *C4  (I)  -  R  A  C  ?  C  (J)  *  S I N  A  (I)  *  (  (FEL  (I,  J+1)  -FEL  (I,  J-1)  ) 

**  (GL  (1+  1  ,  J)  /  S I N  2  A  (1  +  1)  -GL  (1-1  ,J)  /SIN2A  (1-1)  )  )  +  PACP.S  (I)  *C03HA  (J) 

*  *  (  (F  EL  ( I  + 1 ,  J )  -  FEL  (1-1  ,J)  )  *  (GL  (I ,  J  +  1 )  /COS  112  (J  +  i)  -GL  (I,  J-1) 

* /CCS H2  (J-") )  ) 

D2L  =  PGL-GL  ( I , J) 

GL  (I,J)  =GL  (I , J)  +  D2L*?EFG 

IF  (ABS  ( D 1 L)  . LT.  .001. AND. ABS ( D 2 L )  . LT. .001) L  =  L  +  1 
CONTINUE 
GO  TO  3 
DO  23  J=2 , 50 

RFEL=FZL  (I, J  +  1) *C1  (J)  +  FEL (I, J-1)  *C2  ( J) + FEL  ( 1+ 1  ,  J ) *C 3  ( I) 

*+FZL  (T-1  ,  J)  *C4  (I)  -GL  (I  ,J)  *SECH2/C5*  (SIN  11 2  (J)  +COS2A  (I)  ) 

D 1 L=RFEL-?EL  (I , J) 

FEL  (I, J) =FEL  (I, J)  +  D 1 L*  REFF 

RGL=GL  (I  ,  J+  1 )  *C1  (J)  +GL  (I,  J-1 )  *C2  (J)  +GL  (1  +  1  ,  J)  *C3  (I) 

*+GL  (1-1  ,  J)  *C4  (I)  -BACSC  (J)  *SINA  (I)  *  {  (FEL  (I,  J  +  1)  -  FEL  (I,  J  -  1)  ) 

**  (GL  (I,  J)  /3/SINA  (I)  -GL  (1-1  ,  J)  /SIN2A  (1-1  )  )  )  +  RACES  (I)  *C0SHA  (J) 

**  (  (FEL  (1  +  1  ,J)  -FEL  (1-1  ,  J)  )  *  (GL  (I,  J  +  1)  /  JOS  112  (J  +  1)  -GL  (I,  J-1) 

* /COSH 2  (J-1 ) )  ) 

D2L=RGL-GL  (I, J) 

GL  (I ,  J)  =GL  (I ,  J)  +  D2L*R3FG 

IF  (ABS (OIL)  . LT. . 001 . AND. ABS  (D2L)  . LT . . 001 ) L=L+ 1 

CONTINUE 

CONTINUE 

I FLAG =2 

GO  TO  70 


■ 
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?  1 


4  1 


40 
4  2 


4  3 


1  5 
45 

70 

7 

5 

100 


DO  45  J J=2 , 50 
J  =  5 1  +  1  - JJ 
DO  15  1 1- 2 , 3  C 
1=3 1 +1 -II 
IF  (1-30)  40,4  1,40 

?FEL=FEL  (I ,  J  +  1 )  *C  1  (J)  +  ?£L  (I,  J-1)  *C2  (J)  +  FEL  (I  + 1  ,  J ) *  C3  (I) 

*  +  FEL  (1-1  ,  J)  *C4  (I)  -GL  (I,J)  *SECH2/C5*  (SINH2  (J)  +COS2A  (I)  ) 
D1L=3F2L-FSL  (I, J) 

FEL  (I, J) =FEL  (I,  J)  +D1 L*REFF 

RGL=GL  (I,J+1)*C1  ( J) +GL  (I,  J-1)  *C2  (J)  +GL  (1  +  1  ,  J)  *C  3  (I) 

*  +  GL  (1-1  ,  J)  *CU  (I)  -RACSC  (J)  *SINA  (I)  *  (  (FEL  (I,J+1)-FEL(I,J-1)) 

**  (GL  (I,  J)  /B/SINA  (I)  -GL  (1-1,  J)  /SI  N  2  A  (1-1)  )  )  +EACR3  (I)  *COSHA  (J) 

**  (  (FEL (1  +  1 ,  J) -FEL  (1-1 , J) )  * (GL (I, J+1) /COSH 2 (J  +  1  )  -GL (I, J-1) 

* /COS H2  (J-1 ) )  ) 

D2L=RGL-GL  (I, J) 

GL  (I,  J)  =  GL(I,J)+D2L*REFG 

IF  (ABS(DIL)  . LI. . 001. AND. A3S  (D2L)  . LT. . 001 ) L=L+1 
GO  TO  15 
IF  (1-2)  4  2,43,42 

?I  EL=FEL  (I,  J  +  1)  *C1  (J)  +  FEL  (I,  J-1)  *C2  (J)  +  FEL  (1  +  1  ,  J)  *0  3  (I) 

*  +  FEL  (1—1 r  J ) *  C  4  (I)  -GL  ( I ,  J )  *  SEC  H2/C5  *  (SINH2  (J)  +COS2A  (I)  ) 
D1L=RF3L-?3L (I, J) 

FEL  (I, J) =?EL  (I , J) +D1 L*RE?F 

RGL  =  GL  (I,  J  +  1)  ac1  (J)  +  GL  (I,  J-1)  *C2  ( J)  +GL  (I  + 1  ,  J)  *C3  (I) 

*  +  GL  (1-1  ,  J)  *C4  (I)  -EACF.C  (J)  *SINA  (I)  *  (  (FEL  (I,  J+1)  -FEL  (I,  J-1)  ) 

**  (GL  (1+1 , J) /S IN 2 A  (I+1 )  -GL  (1-1  , J)  /SI N 2 A (1-1)  )  )  +RACS3  (I)  *COSHA  (J) 
**  (  (FEL (1  +  1 , J )  -FEL  (I- 1 , J)  )  *  (GL  (I, J  +  1) /COSH 2  ( J  +  1 )  - GL ( I , J -  1 ) 
*/CCSH2 (J-1)  )  ) 

D2L  =  RGL- GL  (I, J) 

GL  (I,  J)  =GL  (I,  J)  +D2L*EE?G 

IF  (A3S  (D1L)  . LT. . 001 . AND. ABS  (D2L)  . LT. .001) L  =  L+1 
GO  TO  15 

?FEL=FEL  (I,  J  +  1)  *C1  (J)  +  FEL  (I,  J-1)  *C2  (J)  +  FEL  (I+1,J)*23(I) 

*  +  FEL  (I-  1  ,  J)  *C4  (I)  -GL  (I  ,  J)  *SECH2/C5*  (SINH2  (J)  +COS2A  (I)  ) 

D 1L=RFEL- FEL  (I, J) 

FEL  (I, J)  =?EL  (I, J)  +D1  L*REF ? 

RGL  =  GL(I,  J  +  1)  *C1  (J)  +GL  (I,  J-1)  *C2  (J)  +GL  (1  +  1  ,  J )  *C3  (I) 

*  +  GL  (1-1  ,  J)  *CU  (I)  -RACEC  (J)  *SINA  (I)  *  (  (F3L  (I,  J  +  1)  -  FEL  (I,  J-1)  ) 

**  (GL  (!+"  ,J)  /SIM  2  A  (1+1)  -GL  (I,  J)  /B/SINA  (I)  )  )  +RACPS  (I)  *COSHA  (J) 

**  (  (FEL (1+1 , J) -FEL  (1-1  , J)  )  *  (GL  (I , J  +  1 ) /COSH 2  (J  +  1 )  -GL (I, J-1) 

* /COSH 2  ( J-1 ) )  ) 

D2L=PGL- GL  (I, J) 

GL  (I,  J)  =GL  (I,  J)  +D2L*EEFG 

IF (ABS (D1L)  .  LT. .001. AND. ABS (D2L)  .LT. .001) L  =  L+1 
CONTINUE 


CONTINUE 
IFL AG=1 
DO  7  I=2,?0 

GTT=  (8.  *FEL  (I,  2)  -  FF.L  (1 , 3)  )  *COSHS/2  .  /A/A/ (S  INHS  +COS2  A  (I)  ) 
GL  (I,  1)  =  (GTI-GL  (1,1))  *RSB+GL  (1,1  ) 

CONTINUE 

IF  (L-1421)  6,  100,  100 
CONTINUE 

WRITE  (6, 500)  RE, AS 
WRITE  (6,501)  N 
DO  9  1=1  ,31 

WRITE  (6, 502)  (FEL  (I, J)  , J= 1 , 1 5) 


■ 


* 


9 


102 


WRITE  (5,999)  (FEL  (I,J)  r  J=  1,51) 

CQNTIN  02 
DO  13  1=1,31 

WRITE  (6,504)  ( r  E  L  ( I , J ) , J  =  1  ,  5 1  ,5) 

I  3  CONTINUE 

WRITE  (6,503) 

DO  11  J  =  1 , 51 
DO  1 1  1=2,30 

VO  (I, J)  =  GL  (I, J) * COS HZ /COS  HA  (J) /SIN A  (I) 

II  CONTINUE 
DO  12  1=1,31 

WRITE  (6, 502)  (VO  (I,J)  ,  J  =  1  ,15) 

WRITE (9 ,999)  (VO (I , J) , J=1 , 51) 

12  CONTINUE 

DO  14  1=1 , 

WRITE  (6,504)  (VO  (I , J)  ,J=1,51,5) 

14  CONTINUE 

WRITE  (6,505)  L 
CDS K 1=0. C 
DO  110  1=1,29,2 

110  CDSKI=SINA  (I)  **2*VO  (1 ,  1)  +  4.  *SINA  (1+  1)  **2*VO  (1+ 1  ,  1)  +SINA  (1  +  2)  **2 
**VO (1+2, 1) +CD5KI 

CDSKI=CDSKI*8. 0*SINHZ/COSHZ/RE/3  . *B 
P1=0.0 

DO  IV  J  =  1  ,49,2 

P1  =  P1+  (  (VO  (2  ,  J)  *u.  0-VO  (3,  J)  )  +4.  *  (VO  (2,  J+ 1 )  *4. -VO  (3,  J+1)  )  + 

* (VG (2,J+2) *4. -VO (5 , J+2)) ) /2.0 

III  CONTINUE 
P1=1.0+3./?S*A/3./B*P1 
P  2 ( 1 ) =0 . 0 

SINA (33) =0.0 
COSA  (33 ) =0.0 
PT  (3  3)  =0.0 
DO  112  1=1 ,2 9,2 
11=1+2 

P2  (IT)  =  (  (-VO  (1, 1)  *3. 0+4.0*VO  (1,2)  -  VO  (I,  3)  )  +  (-VO  (1+1  , 1)  *3. 0 
*  +  4.0“ VO  (1  +  1 ,2)  -  VO  (1+1  ,  3)  )  *4.0  +  (  - VO  ( 1  +  2 ,  1)  *3.0 

*+4.0* VO  (1  +  2,  2)  -  VO  (1+2,  3)  )  )  /A/2.0+  (VO  (I,  1)  +V9  (1  +  1  V)  *4.  O+VO  (1+2,  1)) 
**SINHZ/C05aZ+?2 (I) 

112  ?T(I!)=?1+?2(II)*4.  C/RE/3  .  *D 

PT  (1)  =P1 
C  D  F  0  P.  =  0 . 0 
DO  113  1=1,29,4 

1  1  3  CDFOP=CDFOR+PT  (I)  *COSA  (I)  *SINA  (I)  *2. +  8.*  COSA  (I  +  2  )  *  PI  ( I  +  2) 

**SINA  (1  +  2) +2. *PT  (1  +  4)  *SINA (1+4)  *COSA (1+4) 

CDFOR=CDFOP.*B/1 .  5 
CDTTT=CDSKI  +  CDFOR 
WRITE  (6 , 506) 

WRITE  (6,507)  RE ,  CDSKI ,  C  DFOF. ,  CDTTT  ,  PI 
WRITE  (6,502)  (PT  (I)  ,1=1,31,2) 

AMU=1 . 6  6 7E- 4 
A  IF;LO=.  9267F-3 
CRYLO  =  . Q  2 

RLRSQ=3 . *  A  MU  *  *  2 / A I R  LO/  (CR YLO- AIPLO) *PE**2/3  2 . *CDTIT/980 . /AP 
R  L  F  =  R  L  R  S  Q  *  *  (1  .  /3.  ) 

V  T  E  R  M  =  A  a  U  *  R  E  /  2 .  /?.  L  F.  /A  I  PL  0 
R  L  R=  RL  R  *  1 0  0  0  0 . 


* 
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W 7 1 T E  (5 , 50  3) 

WRITE  (6  ,  5C9)  A  MU  ,  AI  FLO ,  CRY  LO,  A?. 
'WRITE  (5,510)  RLR,  VTFRM 


9  99 

FORMAT  (13E9. a) 

500 

FOPMAT (//, 'REYNOLD  NUMBER  =’,F6.2,',  AXIS  PATIO 

-  1 

9 

F  6 .  3) 

501 

FORMAT (//, 15, '  HERAT ’ON S  REQUIRED') 

5  02 

FOPMAT (/, 1 5F8.  3) 

5  03 

FORMAT (//,'  THE  YORTICITY  FIELD') 

5CU 

FORMAT (/, 11F10.  U) 

5  05 

FORMAT (//,'  TOTAL  GOOD  VALUES  =  ',15) 

506 

FORMAT  (  *  1  '  ,  '  F.EYN  OLD  NO’  ,  9X,  '  CDSKIN'  ,  1  0  X  ,  '  CD  FOR  M  ' 
* ' CDTTT '  , 1 3  X  ,  ’ PO ’ ) 

,10X 

t 

507 

FOPMAT (/, 5E1 3. 5} 

503 

FORMAT (//, 10X, ’VISCOSITY’ ,?X, ’AIR  DENSITY* ,8X, 'C 
*  ?  X , ’ A  XI S  RATIO') 

RYS 

DENSITY ' , 

^09 

FORMAT (/, 4E20. 5) 

5  10 

FORMAT (//,' RADIO S=* , F8 .2,  ’TERM.  VEL=’ , F3. 4) 

STOP 

END 


■  ' 
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C  CALCULATION  OF  FLOW  FIELD  AROUND  A  SPHERE  ** **** ** **** * * ****** * 

DIMENSION  C  3  (61)  ,C4(61),FE(61,81)  ,  G  (6 1  ,  8  1 )  ,  V  (6  1  ,  31 )  ,F  (61,81) 
*,SINFI(61),C0SFI(61)  ,EX?Z  (81)  ,  EXP2Z  (91)  ,  5  INF  12  (61)  , 

* PT  (65)  , P2  (65) 

RE=.  1 
FF=  1.5 
FG  =  0 . 2 

B  =  3. *  3. 1416/180. 

A  =  .  05 

C  5  =  2 . /A /A  +  2. /B/B 
DO  2  J= 1 ,81 
Z=A*FL0AT  (J- 1 ) 

EXPZ  (J)  =  EXP  (Z) 

EXP2Z  (J)  =EXP Z  (J)  *EXPZ  (J) 

2  CONTINUE 

C1=  (2. -A) /2. /A/A/C5 
'C2=  (2.+A) /2.  /A/A/C5 
SINFI  (1  )  =0. 0 
COSFI  (1 )  =1 . 0 
SINFI2 (1 ) =0. 0 
SINFI (6 1 ) =0. 0 
COSFI  (6  1 )  =1 . 0 
SINFI2  (61)  =0. 0 
DO  3  1=2,60 
FI=B*FL0AT  (1-1) 

SINFI  (I)  =S IN  (FI) 

COSFI  (I )  =COS  (FI) 

SI-NFI2  (I)  =  S I  N F I  (I)  * S I N F I  (I) 

COTFI  =  COSFI (I)  /SINFI  (I) 

C3  (I)  =  (2. -B*COTFI)  /2. /B/B/C5  . 

C4  (I)  =  (2. +B*COTFI) /2 . /B/B /C 5 

3  CONTINUE 

F  1  =  1 . 0/EXPZ  (81). 

F3=F1*F1*F1 

F6=F3*F3 

D=2.*(2.*F6-10.-10.*F3+10.*F1) 

DO  6  J= 1 ,81 
DO  6  1=1,61 

FE  (I,  J)  =  (EXP2Z  (J)  /2.  +.  25/EXPZ  (J)  -.75*EXPZ  { J)  )  *  SINFI  2  (I) 

G  (I,  J)  =1  5.  /D/EXPZ  (J).  *SINFI2  (I)  *  (F3*EXPZ  (J)  **3-  1 . 0)  *2.  0 
6  CONTINUE 

DO  4  1=1,61 

FE  (1,81) =EXP2Z(81)  *SINFI2  (I) /2. 

4  CONTINUE 

DO  15  1=2,60 
DO  15  J= 1 , 81 

F  (I,  J)  =3  (I,  J)  /EXP2Z  (J)  /SINFI  2  (I) 

15  CONTINUE 

DO  5  J=1 ,81 

V  (1  ,  J)  =0.0 

V  (61  ,  J)  =0.  0 
F  (1  ,  J)  =0.0 

F  (6  1  ,  J)  =0.0 


. 

. 
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I FL AG=1 
DO  7  N=1 , 1 00 
L=0  • 

IF  (IFLAG.EQ.  2)  GO  TO  20 
DO  9  1=2,60 
DO  9  J  =  2 , 8 0 

RFE=FE (I,J  +  1)  *C1  +  FE  (I, J-1) *C2+FE  (I+1,J) *C3 (I) +FE (1-1,  J) 
**C4  (I) -G  (I, J)  *EXP2Z (J)  /C5 

RG=G  (I,  J  +  1)  *C1+G(I,  J-1)*C2  +  G  (1+1  ,J)  *C3(I)  +  G(I-1  ,  J)  *C4  (I) 
*-RE*EXPZ  (J) *SINFI  (I) /(8.*A*B) *  (  (FE  (I, J+1 ) -FE  (I, J-1)  ) 

**  (F (1  +  1  ,J)  -F (1-1  ,J)  ) 

*-  (FE  (1+1 , J) -FE  (1-1 , J)  )  * 

* (F (I, J+ 1) -F  (I, J-1)  ) )  /C 5 

D1  =  RFE-FE  (I, J) 

D2=RG-3  (I, J) 

IF  (ABS (D1)  .  LT. .001 . AND.ABS  (D2) .  LT. . 001) L  =  L+1 
FE  (I,  J)  =  FS  (I,  J)  +FF*D1 
G  (I, J) =G  (I, J)  +FG*D2 
F  (I,  J)  =G  (I,  J)  /EXP2Z  (J)  /SINFI  (I) 

9  CONTINUE 

IFL AG=2 
GO  TO  60 

20  DO  21  J J=2 , 80 

J=82- J J 
DO  21  11=2,60 


*-RE*EXPZ  (J) * SINFI  (I) /(8.*A*B)  *  ( (FE  (I, J+1) -FE (I, J-1)  ) 
**(F(I  +  1,J)-F(I-1,J)) 

*-  (FE  (1+1 , J) -FE (1-1, J)  )  * 

* (F (I, J+1) -F  (I, J-1)  ) )  /C  5 
D 1  =  RFE- FE  (I, J) 

D2=PG-G  (I, J) 

IF (ABS  ( D 1 )  .LT.  .001, AND. ABS  ( D 2 )  .LT. .001)  L  =  L+1 
FE  (I,  J)  =FS  (I ,  J)  +FF*D1 
G (I, J) =G (I, J) +FG*D2 
F  (I,  J)  =G  (I,  J)  /EXP2Z  (J)  /SINFI  (I) 

21  CONTINUE 


IFLAG=1 


60  IF  (L.  EQ.  5?*'79)  GO  TO  100 
DO  8  1=2,60 

GF.E=  (8.  *FE  (I,  2)  -  FE  (I,  3)  -  7.  *FE  (I,  1)  )  /2.  /A/A 
G  (I,  1)  =  G  (1, 1  )  +  (GRR-G  (1,1)  )  *FG 
F  (1,1)  =G  (I,  1 )  /5INFI2  (I) 


8  CONTINUE 

7  CONTINUE 


8 


7 


100  WRITE  (6 , 500)  N,  L 
DO  10  1=1,61 

WRITE  (7, 999)  (FE  (I  ,  J)  ,  J  =  1 , 81) 
WRITE  (6,501)  (FE  (I,  J)  ,  J=1 , 61 , 5) 


1  00 


10  CONTINUE 

WRITE  (6,50  5) 
DO  1 4  1=2,60 
DO  14  J=1 , 81 


. 
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V  (I,  J)  =  G  (I,  J)  /SINFI  (I)  /EXPZ  (J) 

14  CONTINUE 

DO  13  1=1,61 

WRITE  (6, 501)  (V  (I, J)  , J=1,61  ,5) 

WRITE  (3,999)  ( V  (I ,  J  )  ,  J  =  1 , 8  1 ) 

13  CONTINUE 

CDSKIN=0. 0 
DO  110  1=1,59,2 

1  10  CDSKIN=  SINFI  (i)**2*V(I,1)+4.*SINFI(I+1)**2*V(I  +  1,1)+SINFI(I  +  2) 

***2*V  (1  +  2,1)  +  'DSKIN 
CDSKIN=CDSKIN*4. 0/RE/3 .*B 
P1=0. 0 

DO  111  J=1 ,79,2 

P1=P1+ (  (V  (2,  J)  *4. 0-V  (3  ,  J)  )+4.*(V(2,J+1)  *4. -v  (3, J+1)  )  + 

*  (V  (2  ,  J  +  1 )  *  4.  -  V  (3 ,  J  +  2 )  )  )/2.0 

111  CONTINUE 

PI =1.0+ 9. /3E*A/3. /E*P1 
P 2  (1) =0. 0 
SINFI  (63)  =0.  0 
C05FI  (6  3)  =0.  0 
PT  (63)  =0. 0 
DO  112  1=1,59,2 
11=1+2 

P2 (II)  =  ( (-V (1, 1) *3.  +  4. *V (1,2) -V  (1,3) )  +  (-V  (1+ 1 , 1 )  *3. 

*  +  4 . * V  (1  +  1,2)  -V  (1  +  1,3)  )  *4 .  +  (  -  V  (1  +  2 , 1 )  *3 .  0 

*  +  4.  *  V  (1  +  2,2)  -  V  (1  +  2,3)  )  J/A/2.+  (V  (1,1)  +  V  (1  +  1,1)  *  4 .  +  V  (1+2, 1)  )  +  P2  (I) 

112  PT  (II)  =P1  +  P2  (II)  *4.  /PE/3.  *B 
PT (1) =P1 

CDFOEM=O.Q 
DO  113  1=1,59,4 

1  1  3  CDFORM=C  DFORM  +  PT  (I)  *COSFI  (I)  *  SIN  FI  (I)  *2  .  +8.  *COSFI  (I  +  2  ) 

**SINFI (1  +  2) *PT (1  +  2) +2. *PT  (1  +  4) *SINFI  (1+4)  *COSFI  (1+4) 
CDFORM=CDFCRM*B/1 . 5 
C  DTTT  =  C  DFOP  N  +  CDS  K I N 
WRITE  (6 , 506) 

WRITE  (6 , 507)  RE,  CDS  KIN,  CDFOF.M,  CDTTT,  PI 
A  M U= 1 . 6  67E-4 
A  I R  L  O = .  9267S-3 
CRYLO=1 . 0 

RLRSQ  =  3. *AMU**2/AIRLO/  (CRYLO- AIRLO)  *RE* * 2 /3 2 . *CDTT T /9 8 0 . 
RLR=RLRSQ** (1. /? . ) 

VTERM=AMU*RE/2. /RLR/AIRLO 
RLR=RLR*10000. 

WRITE  (6,508) 

WRITE (6,509) AMU, AIRLO, CRYLO 
WRITE  (6, 510) PLP, VTERM 

500  FORMAT (//, 15, ’  ITERATIONS  REQUIRED',’,  TOTAL  GOOD  VALUES  =«,I5) 

501  FORMAT (/, 2X, 1 3F9. 4) 

505  FORMAT  ('1') 

506  FORMAT  ('  1 ' ,5X, ' REYNOLD  NO ' , 7X , ' CDSKIN • , 8X , ' CDFORM' , 8X, 

* ' CDTTT' , 9X , ' PO ' ) 

507  FORMAT (/, 5E1 5. 5) 

508  FORMAT (//, 1  OX, ' VISCOSITY ' ,7X, ' AIR  DENSITY' ,  6  X , 'DROP  DENSITY') 

509  FORMAT (/, 3E20. 5) 

510  FORMAT (//,' ?ADIUS=' , F8. 2, '  TERM.  VEL.=',F8.3) 

999  FOPMAT  ( 1 3F9. 4) 

STOP 


:  ■ 
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c 

C  TRAJECTORIES  OF  ICE  CRYSTAL  AND  WATER  DROP 

r~% 

L 

C  RUNGE-KUTTA  METHOD  TO  START  INTEGRATION 

DIMENSION  R  P  L  Z  (  4 )  ,RFSZ  (4)  ,  E  RLY  (4)  ,  R  R  S  Y  (4)  ,PVLZ  (4)  ,RVSZ  (4)  , 

*RVLY  (4)  ,  RVSY  (4)  ,  RALZ  (4)  ,  HAS  Z  (4)  ,  RALY  (4)  ,  RASY  (4)  ,W  (61 , 81)  ,0(61,81)  , 
*Z  (61 ,81)  ,Y  (6 1 , 8  1) 

CONS=2.*9 .267/1 .667 
A  R  =  . 05 

ZI  A=  A  LOG  (SORT  ( (1 .  +AR)/(1  ,-AR)  )  ) 

RE  AD  (7 , 1 00)  REYL , R EYS , VTERL , VTE RS , A XI S L , A XISS 
WRITE  (6,110) 

WRITE  (6,1 00)  REYL, REYS, VTERL , VTERS , A  XI SL , A  XI SS 

READ  (7 , 1 C 1 ) DELTAS ,TMAX,Q1  ,Q22,ZZ,YY 

READ  (7, 102) VLZ1 ,ALZ1 , RSZ1 , RSY1 , VSZl  ,ASZ1  ,RSR, RLR 

WRITE  (6, 1C1) DELTAS, TM A X,Q1 , Q2 , Z  Z  ,  Y  Y 

WRITE  (6, 851 )  VLZ 1  ,  AL7.1  ,  P.SZ1  ,  RSY1  ,  VSZl  ,  ASZl  ,  RSR ,  RLR 

DO  701  1=1 ,61 

READ  (1,999)  (W  (I , J)  , J  =  1 , 81 ) 

701  CONTINUE 

DO  702  1=1 , 61 

READ  (2,99  9)  (U  (I , J)  , J=1 , 81 ) 

702  CONTINUE 

DO  703  1=1,31 

READ  (3,999)  (Z  (I, J)  , J=1 , 51 ) 

703  CONTINUE 

DO  704  1=1 , 3 1 

READ  (4,999)  ( Y  (I ,  J)  ,  J=  1 , 51 ) 

704  CONTINUE 

DO  10  1  =  **  ,61 
DO  1C  J  =  1 , 8 1 
W (I, J) =W (I, J) /. 1C6*VSZ1 
10  U  (I, J) =-U  (I, J)/.  1 06*VSY1 

Z  A  D=  0  .  C  0  0 
DO  20  N=1 ,6 
DELT  =  DELTf*S/1  COO.. 

5C  WRITS  (6,852) 

CONFAC=CCNS*RSR*AXISL*VTERL 

RLZ=0 . C 

RLY=C. 0C0 

VL  Z=  VLZ1 

V  L  Y=  0  .  C  0  0 

AL  Z=  ALZ1 

ALY=0. 0C0 

RSZ=R  SZl 

RS  Y=  R  S  Y 1 

VS Z= VSZl 

VSY=0. 000 

A  S  Z=  A  S  Z 1 

AS Y=0 . 000 

RRLZ  (1 ) =  PLZ 

RRLY  (1) =  RLY 

R VLZ  (1 ) = VLZ 

RVLY  (1) =  VL Y 

RALZ  (1 )  =  ALZ 

RALY  (1)  =  A L Y 

PRSZ  (1) =RSZ 


•  ; 


■ 


. 


FRSY  ( 1 ) ~ B S Y 
RVSZ  (1) =  VSZ 
R  VS  Y  (*’)  =  VSY 
RASZ  (1)  =ASZ 
B  A  3Y  (1 )  =  A  S Y 
DO  1  1=1,3 
T=DELT*  (FLOAT  (I)  ) 

ZFA-ZZ-BLZ+RSZ 
YF  A  =  Y Y-PLY  +  PSY 
ZSA=ZFA/RSR*RLR 
YS  A= Y  FA/PSR  *RLR 

AK1=DELT*Q1*  ( 1  •  -  (VLZ-UDR  (ZSA,YSA,W)  )  ) 

AL1  =DELT*  (-Q1* (VLY-UDR (ZSA, YSA,U) ) ) 

AN1=DELT*  (-Q2* (VSY-UCR (ZFA, YFA ,  Y)  )  ) 

AM  1  =  BELT*  (Q * - Q2* (VSZ-UCR (ZFA, YFA, Z)  )  ) 
ZFl=ZFA+DELT/2. * ( VSZ-VLZ) 

YFl=YFA+DELT/2. * ( VSY- VLY) 

ZS1=ZF1/F.SR*RLR 

YS1=YF1/ESR*PLR 

Q2=Q22*REYS/ (VSZ+AM1/2.-UCR  (ZFA , Y F A , Z ) ) /COMFAC 
AK  2  =  DELT*Q1 *  (1 . - ( VLZ  +  AK1/2. -UDR (ZSl,YSi,W))  ) 
AL2~  DELT* (-Q^  * (VLY+AL1/2.-UDR (ZS1  ,  YSl  ,  U )  )  ) 

AM2  =  DELT*  (Q1-Q2*  (VSZ  +  AM1/2.  -UCR  (ZF1  ,  YF1  ,  Z)  )  ) 

AN 2=  DELT* (-Q2 * ( VS Y  + AN1/2 . -UCR (ZF1 , YF1  , Y  )  )  ) 
ZF2=ZFA  +  DELT/2. * (VSZ-VLZ+ (AHl-AKl ) /2. ) 

YF2= Y  FA  +  DELT/2 . * (VSY-VLY+ (AN1-AL1 )  /2.  ) 

ZS2=ZF2/KSR*PLR 

YS  2= YF2/RSR*FLR 

Q2=Q22*REYS/ ( VSZ+ AM2/2. -UCR (ZF1 , YF1 ,Z) ) /CONFAC 
AK3=DELT*Q1  *  (1  .  -  (VLZ  +  AK2/2.  -UDF.  (ZS2,  YS2  ,  W)  )  ) 

AL  3  =  DELT*  (-Q* * ( VLY  +  AL2/2 . -U DP (ZS2 , YS2 , U) ) ) 

AM 3=  DELT*  (Q1-Q2* ( VSZ  +  AM2/2. -UCR (2F2,YF2,  Z)  )  ) 

AN  3= DELT*  (-Q2* (VSY  +  AN2/2 . -UCP  (ZF2 , YF2 ,Y) )  ) 
ZF3=ZFA+DELT*  (VSZ-VLZ+ (AM2- AK2) /2 . ) 

YF 3= YFA* DELT*  (VSY-VLY+ (AN2-AL2) /2.) 

ZS 3^ZF3/FSR*RLR 
YS3=YF3/RSR*RLR 

Q2-Q2  2*R  EYS/ ( VS Z  +  A M3) -UCR (Z F2 , YF2 , Z) ) /CONFAC 
AK4=DELT*Q1*  (1 .  -  (.VLZ  +  AK2-UD  R  (ZS3  ,  YS3  ,  W)  )  ) 

AL4  =  DEL?* (-Q1 * ( VLY+AL3-UDR (ZS3,YS3,U)  ) ) 

AM  4=  DELT* (Q1-Q2* (VSZ+AM3-UCR(ZF3,YF3,Z)  )  ) 

AN  4=  DELT*  (-Q2*  (VSY  +  AN3-UCR (2F3,YF3,Y)  )  ) 
RLZ-RLZ+DELT* (VLZ+ (AK1 +AK2+AK3) /6. ) 

PL Y=R  LY  +  DELT* ( VLY  + (AL1+AL2+AL3) /6 . ) 
RSZ=RSZ+DELT*  (VSZ+(AMl+AM2  +  AM3)/6.) 
RSY=RSY+DELT* (VSY+ ( A  N 1 +  A  N  2  +  A  N  3 ) /6 . ) 

V  L  Z=  VLZ+ (AK1+2.  *AK2+2.*AK3+ AK4) /6. 

VLY=VLY+ (AL1 +2. *AL2+2.*AL3+ AL4) /6 . 

VSZ=VSZ+  (AMI  +2.  * A  M2  +  2.  * A M3  +  AM4 ) /6 . 

VS  Y=  VS  Y+  (AN1  +2.  *A  N2 +2  .  *  A  N3  +  AN4)  /6. 

DP=YY+RS Y-RLY 
DH=ZZ+RSZ-RLZ 
DHS=  DH/RSR*RLR 
DPS=D  P/RSP*R  LR 

Q2=Q22*PEYS/  (VSZ-UCR  (ZF3,YF3,Z) ) /CONFAC 
ALZ  =  Q 1  *  (1.-  (VLZ-UDR  (  DHS,  DPS  ,W)  )  ) 

AL Y=- Q1 *  (VLY- UDR ( DHS , DPS,  U)  ) 

ASZ=Q1-Q2*  (VSZ-UCR (DH,DP,Z)  ) 

A S Y=-Q2* (VSY-UCR (DH, DP, Y)  ) 

WRITE  (6,888)  A K 1 ,AK2,AK3,AK4 ,AL1,AL2/AL3,AL4, 
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SAM1 , A M2, AM  3, AM4 , AN1 , A N2 , AN3 , AN 4 
8  88  FORMAT  (4 F20. 10) 

WRITE  (6,853)  DP, DH , RLY , RLZ, RSY, RSZ ,T 
WRITE  (6,851)  VL2,VLY,  VSZ,  VSY  ,ALZ,.ALY,  ASZ,  ASY 
RRLZ  (1  +  1) =  RLZ 
RVLZ  (1  +  1 ) =  V  L  Z 
RALZ  (I  +  1)  =  A  LZ 
PRLY  (I  +  "1 )  =RLY 
RVLY  (1  +  1) =VLY 
RALY  (I  +  1 ) = ALY 
F RSZ  (1  +  1)  =r.sz 
R  VSZ  (I  +  1 )  =VSZ 
R ASZ  (1  +  1) “ASZ 
RRSY  (I  +  1 ) =RS Y 
RVSY  (1  +  1) =  V  S  Y 
RASY  (1  +  1 )  =  A S Y 
1  CONTINUE 

WRITE  (6,852) 

C  HAMMING'S  PREDICTOR-CORRECTOR- MODIFIER  METHOD 

T=4. *DELT 
TK-4 . *DELT 
VLZ3P=PVLZ  (4) 

V  L  Z3C  =  RVLZ  (4) 

VLY3P=RVLY  (4) 

VL Y3C  =  RVLY  (4) 

VSZ3P  =  PVSZ  (4) 

VSZ3C=  R VS Z  (4) 

VS Y3  P-RVS Y  (4) 

VS Y3C  =  RVSY  (4) 

3  VLZ4P=RVLZ  (1)  +4./3.*  (2.*RALZ  (2)  -  FALZ  (3)  +2.*RALZ  (4)  )  *DELT 

VS  Z4  P  =  ?VSZ  (1)  +4./3. * (2. *RASZ (2) -RASZ  (3)  +2.*RASZ  (4)  ) *DELT 
VLZ4M  =  VLZ4P-1 12./121 . *  (VLZ3 P-VLZ3C) 

VSZ4M=VSZ4°-1 12./121 . *  (VSZ3P-VSZ3C) 

VLY4P  =  RVLY  (1)  +4./3.*  (2.* RALY  (2) -RALY (3)  +2.* RALY  (4) )  * D ELT 
VSY4P=RV5Y  (1)  +4./3.*(2.*RASY(2)-RASY(3)  +2.*RASY(4))*DELT 
VLY4M=VLY4P-1 12./121 . * (VLY3  P-VLY3C) 

VSY4M=VSY4P-1 12./121 .* (VSY^P-VSY3C) 

RLZ4  =  .  ‘*25*  (9 . *R  RLZ  (4)  -RPLZ(2)  +3.*DELT*  (VLZ4M  +  2.  *RVLZ  (4)  -  RVLZ  (3)  )  ) 

RSZ4=.125*(9.*RRSZ(4)  -PPSZ  (2)  +3. * DELT * (VSZ4M+2. *RVSZ (4) -RVSZ  (?) )  ) 

RLY4  =  . 12  5*  (9 . *RRL Y (4) -RRLY (2)  +3 . *DELT* (VLY4M+2. *RVLY (4) - RVLY  (3) ) ) 

RSY4  =  .125*(9.*RRSY(4)-RRSY(2)*3.*DELT*(VSY4?!+2.*RVSY(4)-RVSY(3)j) 

DH=ZZ+RSZ4-RLZ4 

DP=Y Y+RSY4-RLY4 

DHS=DH/RSR*RLR 

DPS=DP/RSF.*RLR 

Q2=Q22*REYS/ ( VSZ4M-UCP (DH,DP, Z) ) /CONE AC 
ALZ4=Ql* (1 . - (VLZ4 M-UDR (DHS, DPS, W)  )  ) 

ALY4  =  -Q1*  (VLY4M-UDR (DHS,DPS,U)  ) 

ASZ4  =Q1 -Q2*  (VSZ4M-UCR  (DH,DP,Z)  ) 

ASY4  =  -Q2* (VSY4M-UCR (DH,DP,Y)  ) 

VLZ4C=.  12  5*  (9  .*RVLZ  (4)  -  RVLZ  (2)  +  3.  *DELT*  (ALZ4  +  2.  *R  A  LZ  (4)  -  RALZ  (3)  )  ) 
VSZ4C=.125*(9.*RVSZ(4)-RVSZ(2)+3.  *  DELT*  (ASZ4+2. *R ASZ ( 4) - RASZ  (3))  ) 
VLZ4F=VLZ4C+9./121 . * (VLZ4P- VLZ4C) 

VSZ4F=VSZ4C+9./121 .* (VSZ4P- VSZ4C) 

VLY4C=. 125* (9.* RVLY (4) -RVLY (2)  +3 . *DELT*  (ALY4  +  2. *RALY (4) - RALY  (3)  )  ) 
VSY4C=.125* (9  .  * R  VS Y (4)  -PVSY  (2)  +3. *  D  ELT*  (ASY 4 +  2. *RASY(4)-RASY(3))) 
VLY4F  =  VLY4C  +  9./1 21 . *  ( VL Y 4 P- VL Y 4 C) 

VSY4F=VSY4C  +  9./121 .*  ( VSY4P- VSY4C) 

RLZ4  =  . 12  5* (9.*RRLZ  (4) -RRLZ (2)  +  3 . * DELT * ( VLZ4 F+ 2 . *FVLZ  ( 4) - RVLZ  (?) )  ) 
RSZ4  =  .125*  (  9 . *R RSZ  (4)  -  RRSZ(2)  +  3 . *  DE  LT  *  (VSZ4F+2.  *  PVSZ  (4)  -RVSZ  (3)  )  ) 


■ 


RLY4=.125*(9.*RRLY  ( 4 )  -RRLY  (2)  +  3  .  *  D  E  LT  *  (VLY4F+2.*RVLY  ( 4 )  -RVLY  (?)  )  ) 
RSY4=.125*  ( 9 . *F  P  S  Y  (4)  -  R  FS  Y  (2)  +3.*DELT*(V5Y4F+2.*PVSY  ( 4 )  - RV  SY  (?)  )  ) 
DP=YY+RSY4-RLY4 
DR- DP-  RLR 

IF  (DR.GT.O.O) GO  TO  510 
DH-ZZ+RSZ4-RLZ4 
DIIS=DH/RSR*RLR 
DPS=  DP/RS  R*  R  LR 

Q2=Q22*REYS/(  V5Z4F-UCP.  (DH,DP,Z)  )  /CONFAC 
ALZ4=Q1* (1 . - ( VLZ4 F-UDR (DHS, DPS, K) ) ) 

ALY4=-Q1* (VLY4F-UDR (DHS, DPS ,0) ) 

ASZ4=Q1-Q2*  ( VSZ4F-UCP .  (DH,  DP,Z)  ) 

ASY4=-Q2* (VS Y4F-UCR  ( DH , DP , Y) ) 

DPOV=DP/C OSH  (ZI  A) 

IF  (DPOV . GT.  1 . )  DPQV-1 . 

ZAD=COS  (ARSIN  (DPOV)  )  *SINH  (ZIA)  +FLOAT  (2*N-1)  *RSR 
IF  (DH-ZAD)  50 1 ,501  ,  51  0 

510  IF  (TMAX.LE.T) GO  TO  501 
IF  (DH.GT. 1 .0) GO  TO  521 
IF  (TK-.'O)  521 ,522, 522 

522  WRITE  (6,850)  DP, DH , RLY4 , PLZ4 , RSY4 , FSZ4 , T 

WRITE  (6, 851)  VLZ4F, VLY4F, VSZ4F, VSY4F, ALZ4 , ALY4 ,ASZ4 , ASY4 
TK  DELT 

521  T-T+DELT 

TK-TK+DFLT 
DO  4  1-1,3 
RALZ  (I) =RALZ  (1+1) 

RVLZ  (I) =  FVLZ  (1+1) 

FRLZ  (I) = R R L Z  (1+ 1 ) 

RALY  (I)  =  P A L Y  (1+1) 

RVLY  (I) =FVLY  (1+1) 

RRLY  (I) =RRLY  (1+1) 

RASZ  ( I )  =  R A S Z  (1+1) 

RVSZ  (I)  -RVSZ  (1  +  1) 

RRSZ  (I)  =  FRSZ  (1+  1 ) 

RASY  (I) =P ASY  (1+1) 

RVSY  (I) =F VS Y  (1  +  1 ) 

FRSY  (I) =RRSY  (1+1) 

4  CONTI NUF 

RALZ  (4) -ALZ4 
RVLZ  (4)  =  VLZ4F 
VLZ3  P= VLZ  4P 
VLZ3C= VLZ4C 
RRLZ  (4)  =RLZ 4 
RALY  (4) =ALY4 
RVLY  (4) -VLY4F 
VL  Y3 P=VLY4P 
VL Y3C-VLY4C 
RRLY  (4 )  =FLY 4 
RASZ  (4) =ASZ4 
RVSZ  (4)  =  V  SZ4  F 
VSZ3P=VSZ4P 
VSZ3C= VSZ4C 
RRSZ  (4) =PSZ4 
RASY  (4) =ASY4 
RVSY  (4) -VSY4F 
VSY3P=VSY4P 
VS Y3C-VS Y4C 
RRSY  (4) =RSY4 
GO  TO  3 


,  ■  ... 
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501  YAD=Y  Y  +  PSY4-RLYU 

WRITE  (6,600)  ...ja 

WRITE  (6,7C0) YY, YAD,ZAD,RSY4 ,RLY4,FSZ4 ,RLZ4 
WRITE  (6,650)  T 
YY=Y Y-RSR/2. 

20  CONTINUE 

100  FO  RM  AT  (6F1 1 . 3) 

110  FORMAT  (2Xr*  RE  CRYSTAL  RE  DROPLET  V  CRYSTAL  V  DROPLET  R 
*  n  DROPLET’) 

101  FORMAT  (2F10. 1 ,2 El 5.4,2F10.2) 

102  FORMAT  (8F5.  4) 

6  00  FORMAT  (//,6X, ’RSY1 ’  ,8X,  ' YAD'  ,8X, ' Z AD'  , 8 X , « RS Y4 '  , 8 X , ’ RLY4 • , 

$  3 X , ' RSZ4 ’  ,8X  ,  ’ RLZ4  ' ) 

999  FORMAT  (1  3F9. 4) 

850  FORMAT  (7F15. 5) 

851  FORM  AT  (8F1 5. 9) 

852  FORMAT  (//) 

7  00  FORMAT  (/,  7  FI  2. 4) 

650  FORMAT  (/,'  TIME  USED  =  ',F8.3,'  SECONDS') 

STOP 

END 

FUNCTION  UCR  (RZ ,R Y , P) 

DIMENSION  P  (31, 51 )  , B (4) 

AR=. 05 

A 1  =.  1 

B1=6. *3. 1416/180. 

Z  I  ANALOG  (SQRT  (  (1  .  +AR)  /  (I  ,-AR)  )  ) 

A  Z  =  A  B  S  ( P  Z ) 

A Y-ABS  (RY) 

IF  (A  Y.GT.  0.0)  GO  TO  1 
ZIN=ALOG  (AZ  +  SQRT  ( AZ*AZ  +  1 . ) ) 

A  N  G=  0 . 0 
GO  TO  49 

1  IF  (AZ. GT. 0. 0) GO  TO  2 
ZIN=0. 0 

ANG= 1 . 5708 
GO  TO  49 

2  Y2=AY*AY 
Z2=AZ*AZ 
C=Y2+Z2-1 . 

B2P4AC=C**2  +  4.  *Z2 

X=  (C+SyFT  (B2P4AC)  ) /2 . 

SI NH  AN=SQRT (X) 

IF  (SINHAN-0. )  40,40,41 
40  ZI N=0 . 0 

A  N  G=  0 . 0 
GO  TO  49 

4  1  SRSQP1=SQBT  (SINHAN*SINHAN  +  1 .) 

ZI N= ALOG  (SINHAN  +  SRSQP1) 

COSAN=AZ/SINH AN 
IF  (COSAN. GE. 1 .0)  GO  TO  56 
ANG= ARCOS  (COSAN) 

GO  TO  49 
56  ANG=0. 0 

4  9  IF  (RZ.GT.C. 0)  GO  TO  50 

A  NG=  3 .  14159 6- A NG 
50  ZIN2  -  ( ZI N- Z I A ) /Al +1 . 

Z I N2 =  (ZIN-ZIA) /Al +1 . 

ANG2= AMG/B1 + 1 . 

NZIN=ZIN2 
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112 
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114 

115 

1  20 

121 

122 

1  23 

C 

130 

C 

1  40 

145 

3 

1 

2 

4  9 
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NANG-AFG2 
IF  (NZIN-  1)  1 1  4 , 1 


0,11 


IF  (NZIN- 50)  1*2,120, 115 
IF  (NANG-1)  121,130,113 
IF  (NANG-30)  1  40,1  30,1  23 
N  Z I N  =  1 
GO  TO  12C 
NZ IN=  50 


IF  (MANG-1)  121,130,122 
N ANG= 1 
GO  TO  130 

IF  (NANG-31)  1  30,1  23,1  23 

NA  NG  =  30 


BI-LINEAP  INTER POL AT ION- EXTRAPOLATION 

NZ IN  Pi =NZI N+ 1 

N ANGP1=N  ANG  +  1 

DZIN  =  ZIN2-FLOAT  (NZIN) 

DA NG~ ANG2-FLO AT  (NANG) 

UCR=P  (NANG,  NZIN)  +  (P  (NANG,  NZIN  Pi)  -P  (NANG,  NZIN)  )  *DZ  IN 
$+  ( P (NANG PI , NZIN)  -P  (NANG, NZIN)  )  *DANG+ ( P (NANG , NZIN) 

$  +  P (NANGP1 ,NZINP1)  -P  (NANGP1 , NZIN) -P  (NANG , NZINPl )  ) *DANG*DZIN 
GO  TO  3 

3 1 - QU A  DP  AT I C  INTERPOLATION 
D  Z I N  =  Z I N  2 - F  L  O  AT (NZIN) 

D  A  N  G~  A  NG  2-  FLO  AT  (NANG) 

NZ INP2=N ZI N+ 2 
NZINPl =NZ IN  +  1 
NZINM1=NZIN-1 
FCT=  (DZIN**2-DZIN)  /4. 

FET= (DAN D**2- DANG) /4. 

DO  145  L  = 1 , 4 
N  =  NA  NG-2  +  L 

3  (L)  =P  (N,  NZIN)  +  (P  (N,  NZINPl )  -P  (N,  NZIN)  )  *DZ  IN  +  (P(N,  NZIN  Ml) 

$  +  P (N , NZINP2) -P (N, NZIN) -P  (N, NZINPl )  ) *FCT 
CONTINUE 

UCR=B  (2)  •*-  (3  (3)-  B  (2)  )  *DANG+  (B  (1 )  +B  (4)  -B  (2)  -B  (3)  )  *FET 

RETURN 

END 

FUNCTION  U  DP  (RZ,RY,P) 

DIMENSION  P  (61 ,81 )  ,3B  (4) 

Bl=3./180.*3. 1416 
A1=. 05 
AZ^ABS  (RZ) 

AY  =  ABS  (P  Y) 

IF  (AY.GT.0. 0)  GO  TO  1 
Z S  P=  ALOG  (AZ) 

A  N  G=  0 . 0 

IF  (AZ. GT. 0. 0)  GO  TO  2 
ZSP=0. 0 
AN  G=1 . 57C  8 
GO  TO  49 

ZS  P=  A  LOG  (SQRT  ( AZ* AZ+ AY* AY) ) 

IF (ZSP.LF. 0. 0) ZSP=C. 0 
ANG=  3. 1 4 1 6- A TAN (AY/AZ) 

IF  (RZ. LT. 0.0) A  NG=  3 . 1 41 6-ANG 
ZIN2=ZSP/A1 +1 . 

ANG2= ANG/B1 + 1 . 

NZIN-ZIN2 
NA  NG  =  A  NG2 

IF  (N Z IN- 1 )  114,120,111 


=*fc 


■ 


•  - 1  \  ■ 


4- 


115 


111  IF  (NZIN-80)  1  1  2,  120, 1  1  5 

112  IF  (NANG-1 )  121 , 1  30 , 1  1  ? 

113  IF  (NANG-bO)  140,  130,123 

114  NZIN=1 
GO  TO  120 

115  UD R=  0  .  C 
GO  TO  3 

120  IF  (NANG-1) 121 ,130 ,122 

121  N  A  NG  =  1 
GO  TO  130 

122  IF  (NANG-61)  130,123,  123 

123  N  ANG  =  60 

C  BI-LINEAB  INTERPOLATION-EXTRAPOLATION 

130  NZINP1=NZIN+1 

N  A  N  G  P 1  =  N  A  N  G  + 1 
DZIN=ZIN 2-FLOAT (NZIN) 

DANG=ANG2-FLOAT (NANG) 

UD  R=P  (NANG, NZIN) +  (P  (  N A NG , NZ I N Pi )  -P  (NANG, NZIN)  ) *DZIN 
$+ (P  (NANGP1 , NZIN) -P  (NANG, NZIN)  )  *DANG+ ( P (NANG, NZIN) 

$  +  P  ( N A N G P n , NZINP1)  -P  (NANGP1 , NZIN) -P  (NANG , NZINP1 )  ) *DANG*DZIN 
GO  TO  3 

C  BI-QUADPATIC  INTERPOLATION 

1 40  DZIN=ZIN2- FLOAT (NZIN) 

D A  NG= ANG2-FLO  AT (NANG) 

NZINP2=NZIN+2 
NZ IN  P 1  =  NZIN  +  1 
NZ IN  Ml =NZI N- 1 
FCT- ( DZIN**2 - DZIN) /4 . 

FET=  (DAND**2-DANG)  /4  . 

DO  145  L=1 , 4 
N  =  N  A  N  G  -  2  +  L 

BB  (L)  =P  (N  ,  NZIN)  +  (  P  (N,  NZIN  Pi  )  -P  (N,  NZIN)  )  *DZIN+  (P  (N,  NZIN  Ml ) 

$  +  P  (N,NZINP2) -P (N, NZIN) -P (N, NZINP1)  )  *FCT 
145  CONTINUE 

UDR  =  B  B  (2)  +  (BB  (3)  -EB  (2)  )  *DANG+  (BB  (1)  +BB  (4)  -B3  (2)  -BB  ( 3)  )  *FET 
3  RETURN 

END 


ILE 


* 


